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Studies  of  dynamic  systems  have  shown  that  oscillations 
in  time  and  space  are  related,  both  being  generated  by  non- 
linear, pulsing  behavior  that  is  derived  from  the  mathema- 
tics of  energy  processing.   Similar  mathematics  exist  in 
chaos  theory,  bifurcation  theory,  and  catastrophe  theory. 
Production-consumption  models  that  simulate  pulsing  proper- 
ties of  ecological  systems  are  of  this  class.   This  dis- 
sertation examines  the  spatial  patterns  and  energetics  of 
autocatalytic  and  pulsing  models  as  a  paradigm  for  ecolog- 
ical and  general  systems.   Configurations  were  tested  with 
steady  or  varying  resource  availability  for  ability  of  the 
model  systems  to  maximize  power  as  the  criterion  for  utility 
and  success.   The  spatial  distribution  of  gaps  generated  by 
simulations  was  compared  to  that  observed  in  rain  forests. 

Models  studied  included  (a)  aggregated,  single- 
compartment  autocatalytic  designs;  (b)  parallel  production- 
consumption  design;  (c)  production-consuraption-recycle 


designs;  and  (d)  multiple  cell  spatial  models  each  with  a 
unit  model  but  interconnected  in  different  ways. 

Models  with  autocatalytic  feedbacks  utilized  more 
power  than  the  same  models  with  only  linear  pathways.   Per- 
cent power  used  increased  with  increasing  available  power. 

Production-consumption  models  show  multiple  steady 
states  with  pulsing  behavior  as  a  transition  between  two 
steady  states.   Localized  maxima  of  power  use  occur  during 
pulsing  but  the  overall  power  use  is  related  to  input  power. 

Spatial  patterns  of  production  and  consumption  in 
spatial  models  were  related  to  input  energy  patterns,  the 
degree  of  connectivity  between  the  individual  cells  in  the 
model,  and  the  hierarchical  level  of  intercell  connections. 
Large  variations  in  patterns  were  accompanied  with  small 
changes  in  power  utilized. 

Edges  of  a  spatial  system  can  act  as  a  source  or  sink 
for  energy  depending  on  the  relationship  between  available 
energy  inside  and  outside  the  boundaries  and  the  degree  of 
connectivity  along  the  edges. 

Basic  autocatalytic  production-consuraption-recycle 
models  with  different  spatial  conditions  organize  different 
spatial  patterns  while  generating  near  total  utilization  of 
available  power.   The  wide  variety  of  spatial  patterns 
results  from  dynamic  adaptations  for  maximizing  power  for 
different  spatial  conditions.   The  simulation  results 
resemble  patterns  in  nature  often  attributed  to  random 
indetermi nancy. 
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CHAPTER  1 
INTRODUCTION 

Ecosystems  develop  patterns  in  time  and  space.   Some  of 
these  patterns  are  generated  by  pulsing  oscillatory  proces- 
ses.  What  sorts  of  interactions,  organization  and  structure 
in  an  ecosystem  lead  to  pulsing  behavior,  and  how  does  this 
behavior  affect  the  use  of  energy?   What  types  of  spatial 
patterns  develop  when  ecosystems  are  influenced  by  pulsing 
in  time  and  space?  What  are  the  energy  implications  of 
different  pattern  forming  processes  in  ecosystems?   What  are 
the  effects  of  pulsing  on  succession,  competition,  frequency 
response  of  producers  and  consumers,  and  coupling  with  ex- 
ternal pulses? 

This  dissertation  uses  general  systems  models  to  ana- 
lyze the  effects  of  pulsing  on  pattern  formation  and  overall 
power  use  as  systems  develop,  build  structure  and  organize 
in  time  and  space.   Simulation  models  using  general  systems 
principles  and  based  on  real  ecosystems  were  used  to  test 
the  role  of  pulsing  behavior  of  consumers  in  organizing  eco- 
systems over  time  and  space.   Data  from  a  tropical  ecosystem 
were  used  to  calibrate  pulsing  and  spatial  models. 


Historical  Perspective 

Previous  Models  of  Pulsing  Patterns  in  Time  and  Space 

In  many  fields  from  chemistry,  physics,  and  biology  to 
astronomy,  there  are  a  variety  of  models,  methods  and  tech- 
niques to  describe  and  study  systems  that  have  discon- 
tinuities or  other  rapid  fluctuations  in  their  behavior. 
Some  of  these  are  catastrophe  theory  (Thom  1975),  bifurca- 
tion theory,  synergetics  (Haken  1977a, 1977b  ,1979) ,  dynamical 
system  theory  (Rosen  1970),  chaos  and  order  (Prigogine 
1980,1984,  and  Schaffer  and  Kot  1985),  pulsing  (LotJca  1920 
and  Odum  1982),  pattern  recognition,  and  morphogenesis 
(Meinhardt  1982).   In  all  of  these,  processes  being  de- 
scribed are  parts  of  nonlinear  thermodynamically  open  sys- 
tems.  Energy  constraints  on  these  types  of  systems  have  not 
previously  been  well  studied. 

In  the  past,  efforts  to  describe  systems  using  clas- 
sical thermodynamics  centered  on  closed  systems  near  equili- 
brium or  open  systems  near  steady  state.   In  such  systems, 
available  energy  is  small.   These  approaches  using  equili- 
brium thermodynamics  could  not  account  for  the  behavior  of 
many  systems  (Odum  1983,  Prigogine  1984,  Schaffer  and  Kot 
1985) . 

Data  with  statistical  anomalies  are  often  difficult  to 
analyze  and  methods  are  sometimes  used  to  minimize  fluc- 
tuations (Piatt  and  Denman  1975).   Systems  that  have  aperi- 
odic behavior,  a  great  deal  of  noise,  or  time  dependent 
changes  in  variance  are  not  well  suited  to  the  normal 


statistical  methods.   These  'unusual  events'  can  be  impor- 
tant in  understanding  how  a  system  works  (Weatherhead  1986). 
Frequency  analysis  has  been  used  for  some  time  to  study 
periodic  behavior  of  systems  (Piatt  and  Denman  1975,  and 
Emanuel,  West  and  Shugart  1978).   Fourier  transformations 
decompose  the  output  or  behavior  of  a  system  into  an  addi- 
tive series  of  sinusoidal  processes.   The  variance  is  parti- 
tioned into  a  set  of  frequencies  that  when  combined  gives 
the  output  being  measured.   Aperiodic  behavior  or  systems 
with  known  nonlinear  components  may  also  be  studied  with 
these  techniques,  but  the  results  are  often  not  useful. 
Some  nonlinear  systems  with  behavior  described  as  'chaotic' 
have  frequency  domain  variance  as  noisy  as  the  time  domain 
variance  (Abraham  and  Shaw  1984a,  1984b). 

Pattern  Formation 

Patterns  in  natural  systems  range  from  the  smallest 
molecular  patterns  of  motion  to  the  placement  of   the  stars 
and  galaxies  in  the  universe.   One  of  the  most  intriguing 
aspects  of  pattern  formation  is  the  similarity  of  patterns 
at  differing  time  scales  and  sizes.   From  a  systems  point  of 
view  this  would  lead  one  to  suspect  that  the  processes  are 
similar  at  each  scale. 

Chemically  reacting  systems  give  rise  to  various  types 
of  patterns  (Bray  1921,  Nicolis  and  Prigogine  1969,  Winfree 
1973,  Haken  1977a,  1977b).  The  Belousov-Zhabotinski  reac- 
tion, which  makes  fascinating  patterns,  is  a  simple 


oxidation-reduction  reaction  involving  malonic  acid,  broraate 
and  a  cerium  catalyst  (Winfree  1973).   An  example  of  the 
time  and  spatial  development  of  this  reaction  is  shown  in 
Figure  la. 

Morphological  development  in  biological  systems  has 
been  studied  and  modeled  by  Meinhardt  (1982).   Patterns  form 
when  autocatalytic  growth  in  a  system  is  combined  with 
lateral  inhibition  (negative  spatial  feedbacl^s)  .   Once  auto- 
catalytic activity  starts,  there  must  be  a  longer  range 
negative  feedbacl<  (spatial  inhibition  of  the  spread  of  this 
autocatalysis)  or  the  whole  system  will  pulse  in  a  burst  of 
autocatalytic  consumption.  This  sets  up  spatial  chemical 
gradients  along  which  morphogenesis  is  thought  to  occur 
(Figure  lb) . 

Hilborn  (1979)  experimented  with  predator-prey  models 
based  on  an  aquatic  ecosystem.  Hilborn's  model  had  100 
spatial  cells  arranged  in  a  linear  chain  with  the  ends 
connected   to  form  a  circle.   Both  predators  and  prey  were 
allowed  to  diffuse  across  cell  boundaries.   The  model  was 
simulated  with  initial  conditions  set  so  that  all  cells  had 
prey  but  only  one  cell  had  a  predator.   The  model  (Figure 
2a)  was  allowed  to  iterate  for  1000  time  intervals,  gen- 
erating the  pattern  seen  in  Figure  2b.   Further  experiments 
showed  that  there  was  no  tendency  towards  equilibrium  in 
longer  runs  of  the  model. 

The  spatial  development  of  insect  eyes  and  insect  legs 
has  been  modelled  by  Ransom  (1981)  using  an  autocatalytic 


Figure  1.   Spatial  patterns  based  on  chemical  reaction 
mechanisms . 

(a)   Spatial  patterns  generated  by  Belousev-Zhabotinski 
chemical  reaction  (Prigogine  1980). 


(b)   Spatial  patterns  generated  by  simulation  model 
used  to  describe  morphogenesis  (Meinhard  (1932) . 


Figure  2.   Hilborn's  (1979)  spatial  model. 

(a)  Energy  diagram  of  individual  cell  model 

Equations  for  simulation  model. 

dX(i)=a*X(i)  -b*X(i)*X(i)  -  (c*X  (i  )  *  Y(  i  )  /  (d+X  (  i  )  )  ) 
+h*X(i+l)  +h*X(i-l)  -  2-»;T*X(i) 

dY{i)=-e*Y(i)  -f*Y(i)*Y(i)  +  (g*X  (  i  )  *  Y(  i  )  /  (d  +  X  (i  )  )  ) 
+k*Y(i+l)  +k*Y(i-l)  -  2*k*Y(i) 

where  i  is  the  number  of  the  subsystem  in  a  linear  loop. 


(b)  Simulation  results  of  linear  series  of  unit  models 
showing  level  of  predator  vs  distance  around  loop. 


SPATIAL   LOCATION 


model.   By  allowing  cells  in  the  model  to  divide  and  migrate 
within  given  constraints,  the  model  developed  patterns  simi- 
lar to  those  in  real  insects.   The  model  allowed  simple 
random  cell  division  with  movement  constrained  to  a  hex- 
agonal direction  away  from  the  center  of  the  cell  division. 

Sergin  (1978,  1979,  1980)  studied  the  oscillatory  be- 
havior of  long  terra  climate  variations  using  models  that 
combine  linear  and  nonlinear  interactions  of  the  heat  cap- 
acities of  the  oceans  and  polar  ice  sheets.   The  period  of 
the  climatological  events  in  these  models  is  on  the  order  of 
10,000  to  100,000  years.   The  model  of  global  temperatures 
varies  in  its  behavior  from  steady  state  to  oscillations 
based  on  small  changes  in  areal  coverage  of  continental  ice 
sheets . 

Pattern  formation  based  on  digital,  rule  based  systems 
has  been  used  to  model  biological  systems.   Examples  such  as 
cellular  automata  (Turing  1952  and  Wolfram  1984)  and  a  'game 
of  life'  (Gardner  1970  and  Poundstone  1985)  generate  complex 
spatial  patterns  from  simple  rules.   The  'game  of  life'  is 
generated  on  an  N  x  N  matrix  where 

1.  Every  active  cell  with  two  or  three  neighboring 

cells  survives  to  the  next  generation. 

2.  Each  active  cell  with  four  or  more  neighbors 

'dies'  from  overpopulation.   Every  active  cell 
with  one  or  no  neighbors  'dies'  from  isolation. 

3.  Each  empty  cell  adjacent  to  three  'live'  neighbors 

gives  birth  to  a  new  cell. 
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Figure  3  is  an  example  of  the  patterns  generated  from  a 
simple  five  cell  seed  (R-pentomino)  during  512  iterations. 
This  pattern  stabilizes  (no  more  deaths  and  no  more  births) 
after  1103  iterations,  although  it  is  an  oscillating  steady 
state.   Individual  subsets  of  the  final  stable  pattern 
oscillate. 

The  'game  of  life'  model  has  some  of  the  features  of 
autocatalysis  (or  cooperative  behavior).   Two  or  three  live 
cells  are  required  for  survival  or  birth  of  new  cells.   It 
also  has  the  feature  of  diffusive  inhibition  because  indivi- 
dual cells  that  move  out  from  a  population  center  can  become 
isolated  and  die.   This  rule-based  system  has  no  energy 
constraint  that  governs  development  and  thus  gives  no  energy 
basis  for  pattern  formation. 

The  common  theme  that  runs  through  these  examples  is 
one  of  combined  interactions  of  autocatalytic  growth  with 
some  form  of  inhibition,  diffusion  or  other  mechanism  for 
preventing  the  autocatalytic  growth  from  spreading  too 
rapidly.   A  concept  that  is  sometimes  misunderstood  or  mis- 
interpreted is  that  the  terms  fluctuation  (Prigogine  1980, 
1984)  and  bifurcation  theory  (Pacault  1977)  refer  to  a 
change  in  the  kinetics  of  reacting  components  of  a  system. 
This  change  in  kinetics  gives  rise  to  the  oscillations  or 
pulses  in  the  output. 

The  models  in  this  dissertation  also  use  combinations 
of  autocatalytic  and  diffusion  (linear)  pathways  to  study 


Figure  3.   The  development  of  spatial  patterns  among  cells 
based  on  simple  r-pentamino  initial  condition  (a)  in  'a  game 
of  Life'  simulation. 
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the  possible  mechanisms  and  energy  consequences  of  pattern 
formation  in  ecosystems. 

Concepts  of  Pulsing,  Patterns  and  Power 

Maximum  Power  in  Systems 

Although  in  the  last  century  Podalinsky,  Ostwald  and 
Boltzman  suggested  energy  use  controlled  system  performance 
(Martinez-aiier  1987),  Lotka  (1922)  made  a  more  definitive 
statement.   He  stated  that  evolution  proceeded  in  such  a 
direction  as  to  make  the  total  energy  flux  through  the 
system  a  maximum  compatible  with  the  constraints  on  the 
system.   He  related  this  to  Ostwald's  (1892)  idea  of  all 
possible  energy  transformations,  that  one  takes  place  which 
brings  about  the  maximum  transformation  in  a  given  time. 

A  theory  of  minimum  entropy  generation  was  put  forth  by 
Prigogine  (Prigogine  and  Wiaurae  1946)  that  a  system  evolved 
toward  a  stationary  state  characterized  by  the  minimum 
entropy  production  compatible  with  the  constraints  on  the 
system.   He  has  since  called  this  a  failure  and  probably  a 
special  case  of  systems  near  equilibrium  (Prigogine  1984). 
Prigogine  (1978,  1980,  1982;  Prigogine  and  Stengers  1984) 
now  deals  with  systems  far  from  equilibrium  that  have 
dynamic  and  oscillatory  behavior.   He  has  not  postulated  any 
definite  theory  about  the  energetic  consequences  of  these 
types  of  systems. 

Odum  and  Pinkerton  (1955)  proposed  that  natural  sys- 
tems tend  to  operate  at  that  efficiency  which  produces  a 
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maximum  power  output,  a  general  restatement  of  Lotka's 
original  idea  of  maximum  energy  flow  but  with  an  important 
distinction.   Odum  (1971,  1982,  1983a,  and  1983b)  further 
clarifies  maximum  power  as  useful  power  where  'use'  is 
feedback  of  the  product  of  energy  use  to  amplify  other 
pathways. 

In  describing  cycles  of  life,  death  and  regeneration, 
Calow  (1978)  has  found  that  although  Lotka's  principle 
holds,  there  seem  to  be  no  a  prior  i  grounds  for  placing 
restrictions  on  how  this  use  of  energy  should  be  achieved. 
He  further  stated  that  selection  would  have  shifted  in  the 
course  of  time  from  one  of  maximizing  speed  to  maximizing 
efficiency.   This  is  a  restatement  of  the  strategy  of  eco- 
system development  utilizing  r  and  K  growth  (Odum  1969). 

Jantsch  (1980)  suggests  than  maximum  engagement  in 
matter  (i.e.,  energy  storage)  and  maximum  process  intensity 
(i.e.,  entropy  production)  are  criteria  for  ecosystem 
stability.  Non-equilibrium  structures  thus  come  about  by 
fluctuations  in  the  mechanisms  which  result  in  modifications 
of  the  kinetic  behavior  of  these  structures. 

Design  for  Maximum  Power 

The  important  question  here  is  how  do  systems  build 
structure  in  order  to  maximize  utilization  of  available 
power.   Odum's  theory  (1971  and  1983)  is  that  by  feeding 
back  energy  (derived  from  structure  that  is  being  built) 
reinforcement  occurs  that  increases  efficiencies  and  energy 
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flow  into  the  structure.   Mechanisms  must  develop  that  build 
structure  to  capture  the  most  energy  possible.   These  feed- 
back structures  then  have  a  prior  energy  use  embodied  in 
them  (emergy,  after  embodied  energy,  of  a  structure  has  been 
defined  as  the  total  amount  of  energy  used  in  developing 
these  structures  (Odum  1983  and  1986)).   This  dissertation 
lool^s  at  some  of  the  possible  Icinetic  pathways  that  feed 
back  to  process  energy  and  the  energetics  of  these  pathways. 

Pathway  Configuration 

A  simple  model  demonstrates  several  ways  in  which  use- 
ful power  can  be  increased  (Figure  4,  see  description  of 
symbols  in  Figure  9) .   This  model  is  a  single  storage  with 
autocatalytic  production  drawing  on  a  flow-limited  energy 
source  (an  energy  source  with  constraints  on  the  pathway, 
limiting  the  amount  of  energy  that  can  be  delivered). 

The  efficiency  of  a  pathway  can  be  increased  if  less 
energy  is  fed  back  to  gain  more  energy.   For  a  simple  auto- 
catalytic system  (Figure  4a  and  4b)  this  can  be  done  by 
either  using  less  energy  to  gain  the  same  inflow  (changing 
the  value  of  K2  in  the  model)  or  by  increasing  the  inflow 
for  the  same  feedback  (increasing  Kl  while  concurrently 
decreasing  K3).   Because  there  are  thermodynamic  limits  on 
any  process,  it  may  not  be  possible  to  improve  designs  to 
increase  energy  flows  beyond  thermodynamic  limits. 

The  first  law  of  thermodynamics  requires  the  conserva- 
tion of  energy.   This  implies  the  following  constraint  on 
the  production  process  of  the  model  (Figure  4). 


Figure  4.   Basic  autocatalytic  model  with  Clow-limited 
energy  source. 

(a)  Diagram  with  kinetic  terms 

dQ  =  K1*JS*Q  -  K2*JR*Q  -  K4*Q 
JR  =  JO  /  (I+KO*Q) 

(b)  diagram  with  Clow  terms 


RO  =  :<0*JR*Q 
Rl  =  K1*JR*Q 
R2  =  K2*JR*Q 
R3  =  K3*JR*Q 
R4  =  K4*Q 
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K4*0 
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K0*Jr*Q  +   K2*Jr*Q  =   Kl*Jr*Q  +  K3*Jr*Q        (1) 
Substitute  R  (flow)  terras  as  abbreviations  for  terms 
in  equation  (1 ) : 

R0     +     R2     =     Rl     +     R3  (2) 

Inputs  of  energy  of  any  process  must  equal  the  outputs. 
Efficiency  is  defined  as: 

Efficiency  =  (Output  of  useful  power ) /Inputs 
or  in  terras  of  our  equation: 

Efficiency  =  R1/(R0+R2)  (3) 

where  R3  is  waste  heat  generated  in  the  process  (re- 
quired by  the  second  law  of  thermodynamics).   Because  R3 
cannot  be  zero,  there  is  a  natural  upper  limit  to  the  effi- 
ciency of  any  process. 

Another  method  to  increase  energy  flow  from  a  flow 
limited  source  is  to  have  multiple  pathways  capture  avail- 
able energy,  each  effective  at  a  different  energy  level 
(Figure  5).   Multiple  pathways  (J1,J2,J3)   use  stored  energy 
to  build  structures  to  capture  available  energy.   A  linear, 
donor-controlled  pathway  (Jl)  requires  little  structure  and 
employs  no  feedback  in  order  to  capture  energy,  but  has 
severe  limitations  (its  efficiency  cannot  change)  due  to  the 
dependency  on  the  energy  source.   An  autocatalytic  pathway 
(J2)  feeds  back  embodied  energy  (structure  built  by  the 
system)  to  draw  in  more  energy.   The  quadratic  pathway  (J3) 
is  a  co-operative  phenomenon  in  which  the  structure  of  the 
system  is  interacting  with  itself  to  feed  back  embodied 
energy  to  draw  in  more  power.   A  system  that  develops  such 


Figure  5.   Basic  multiple  path  model.   Three  input 
pathways  represent  different  feedback  regimes: 
linear  (Jl) ,  autocatalytic  (J2) ,  and  quadratic  (J3) 
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higher  order  feedback  pathways  may  exhibit  a  greater  rate  of 
use  of  available  energy. 

This  added  quadratic  pathway  is  available  to  utilize 
any  energy  left  after  the  efficiency  is  raised  to  the  upper 
limit  for  the  autocatalytic  pathway.  This  is  a  mechanism 
that  can  draw  in  energy  that  would  normally  be  unavailable 
to  the  system.   The  quadratic  pathway  may  have  a  high  cost 
to  develop  and  maintain  this  pathway  but  it  enhances  overall 
use  of  that  extra  energy  by  the  whole  system.   This  may  give 
a  competitive  edge  in  some  circumstances  over  systems  with- 
out higher  order  pathways,  particularly  when  available 
energy  may  be  fluctuating.   Available  power  will  be  in- 
creased by  switching  from  one  pathway  to  the  other  depending 
on  the  energy  source.   Some  pathways  are  more  efficient  at 
low  energy  levels  while  others  are  more  efficient  at  high 
energy  levels,  thus  allowing  such  systems  to  efficiently 
utilize  fluctuating  power  sources. 

Pulsing  and  Patterns  in  Ecosystems 

Succession  and  Disturbance 

Any  climax  state  is  eventually  interrupted  by  disturb- 
ances that  generate  patches  in  which  succession  is  re- 
initiated.  The  gaps  in  a  forest  may  be  generated  by  local 
outbreaks  of  consumers  within  the  forest,  tree  mortality,  or 
outside  disturbances  such  as  fires,  hurricanes,  volcanic 
activity,  and  landslides  (Runkle  1985).   The  role  of  the 
landslide  as  a  gap-forming  mechanism  has  been  described  in 
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both  temperate  forests  (Oliver  1981,  and  Veblin  1985)  and 
tropical  forests  (Garwood  1979,  and  Leigh  et  al.  1982). 

Disturbances  (i.e.,  pulses)  to  an  ecosystem  can  be  gen- 
erated from  within  or  can  come  from  outside  the  boundaries 
of  an  ecosystem  and  may  vary  in  frequency  and  amplitude. 
The  ability  of  an  ecosystem  to  utilize  available  resources 
and  adapt  to  these  disturbances  depends  on  the  storages, 
structures  and  interactions  within  an  ecosystem  (Odura  1983). 
Hierarchical  mechanisms  may  develop  that  capture  and  process 
energy  at  various  levels  and  result  in  utilization  of  energy 
over  a  wider  variety  of  input  levels.   Some  mechanisms  of 
interaction  between  parts  of  the  ecosystem  were  studied  in 
this  dissertation  to  understand  how  systems  may  converge 
energy  transformations  and  feedbacli  controls  to  organize  for 
higher  productivity. 

No  unified  theory  of  succession  presented  to  date  can 
be  regarded  as  widely  accepted  (Anderson  1986).   Horn  (1976) 
wrote   'The  sweeping  generalization  that  can  be  safely  made 
about  succession  is  that  it  shows  a  bewildering  variety  of 
patterns.'   Even  the  definitions  of  succession  are  widely 
varying.   In  this  dissertation  succession  is  regarded  as  a 
dynamic  process  in  which  the  composition  of  an  ecosystem 
changes  through  time,  building  structure  and  processing 
energy.   This  process  eventually  stabilizes  in  a  climax  from 
which  there  is  a  regression  or  loss  of  that  structure  due  to 
disease,  fire,  treefall  or  other  events.   Seeding  from  an- 
other ecosystem  or  from  storages  in  the  soil  from  the 
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previous  ecosystem  regenerates  a  facsimile  of  the  original 
ecosystem  through  a  sequence  of  unidirectional  stages  that 
reaches  a  steady  state  system  called  a  climax.   This  climax 
may  be  arrested  at  some  point  and  in  some  cases  succession 
may  cycle  between  several  stages.   This  definition  is 
broader  than  most  but  is  an  attempt  to  describe  the  whole 
process  instead  of  the  more  narrow  'growth-phase'' 
definition . 

Regression  from  a  climax  state  may  occur  in  several 
ways.   In  some  cases  it  comes  about  as  a  pulse  of  con- 
sumption from  within  the  ecosystem  boundaries  such  as  tree- 
falls,  landslides  or  disease  outbreaks.   It  can  also  come 
about  from  disturbances  from  larger  outside  events  such  as 
hurricanes  or  drought.   The  frequency  and  amplitude  of  these 
disturbances  tend  to  be  inversely  correlated:  larger  dis- 
turbances occur  less  frequently  than  smaller  ones.   This 
phenomenon  is  referred  to  as  a  hierarchy  of  disturbances 
(Bennett  and  Chorley  1978).   The  interaction  of  these  dis- 
turbances along  with  the  internal  fluctuations  may  lead  to 
the  'bewildering  variety  of  patterns'  to  which  Horn  refers. 

Edges 

Ecosystems  can  generally  be  broken  up  into  subsystems 
that  have  uniform  characteristics.   These  subsystems  have 
boundaries  where  the  composition  changes  from  one  particular 
type  to  another.   The  development  of  these  edges  may  occur 
where  differing  types  of  energy  interact  with  ecosystem 
components  to  generate  patches  and  zones  of  transition.   The 
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presence  of  many  spatially  distributed  patches  may  be  due  to 
the  production-consumption  pulsing  of  components  in  the 
ecosystem. 

Hierarchies  and  Patches 

The  frequency  of  disturbance  based  on  internal  cycles 
has  been  shown  to  be  from  200-500  years  in  a  variety  of 
ecosystems  (Emanuel,  West  and  Shugart  1978,  Runkle  1985). 
Distribution  of  disturbances  over  time  varies  from  fairly 
constant  low  amplitude  disturbances  to  long-period,  high 
amplitude  disturbances.   The  successional  changes  due  to 
disturbances  may  be  related  to  the  size  and  scale  of  the 
disturbance  (Peet  and  Christensen  1980,  Peet  1981). 

Brokaw  (1982a,  1982b,  and  1985a)  found  a  hierarchical 
distribution  in  gap  sizes  in  a  tropical  rain  forest  at  Barro 
Colorado  Island  (Figure  6a) .   The  area  per  size  class  is 
plotted  vs.  the  size  class  (Figure  6b).   This  relationship 
may  be  important  in  determining  patch  dynamics.   Brown 
(1980)  suggested  that  size  class  distributions  may  be 
related  to  the  emergy  per  size  class  (the  emergy  per  size 
class  is  also  related  to  the  area  per  class)  .   Brolcaw  calcu- 
lated the  turnover  rate  for  the  forest,  based  on  the  gap 
formation,  to  be  from  85  to  128  years  depending  on  the 
minimum  size  of  the  lowest  class  used. 

Models 
The  simulation  models  used  to  study  ecosystem  behavior 
generally  fall  into  two  classes  (Shugart  1984).  One  of  these 


Figura  6.   Size  class  distribution  of  gaps  formed  in 
tropical  forest  at  3arro  Colorado  (3rokaw  1982) . 

(a)  Distribution  of  gaps  by  diameter  of  gap. 

(b)  Distribution  of  gaps  by  area  in  gap. 
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is  based  on  the  nonlinear  'Lotka-Vol terra  equations'  and 
generally  does  not  include  outside  influences.  The  other 
uses  forced  linear  systems  of  differential  equations  and 
does  have  inputs  from  outside  the  system.   Neither  of  these 
methods  typically  contains  any  spatial  considerations  and 
both  deal  with  systems  near  equilibrium.   Systems  near  equi- 
librium tend  to  move  toward  that  equilibrium  and  are  char- 
acterized by  spatial  uniformity  (Prigogine  1984  and  Field 
1985) . 

In  this  study,  open  non-equilibrium  models  are  de- 
veloped that  combine  non-linear  and -oscillatory  interactions 
between  production  and  consumption  with  outside  forcing 
functions  that  provide  resource  controls.   A  pulsing,  hier- 
archical model  of  production  and  consumption  is  used  to 
generalize  about  succession  and  regression.   Spatial  inter- 
actions generated  by  this  model  are  studied  to  understand 
the  energetic  and  kinetic  basis  for  pattern  formation  in 
ecosytems. 

Gap  Models  and  Patch  Dynamics 

Several  previous  studies  based  ecosystem  models  on 
disturbance  gaps.   The  JABOWA  forest  simulator  model  by 
Botkin,  Janak  and  Wallis  (1972)  keeps  track  of  the  birth, 
growth,  and  death  of  a  group  of  trees  from  seedlings  on  to 
maturity  within  a  certain  gap  size.   Subroutines  are  used 
for  crowding,  shading,  and  response  to  individual  nutrients 
and  energy  sources.   The  simulation  then  allows  the  gap  to 
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develop  a  distribution  of  trees  based  on  all  of  the  input 
parameters.   These  gap  models  generally  do  not   account  for 
any  outside  disturbances  that  generate  gaps. 

Various  gap  models  (Phipps  1979,  Shugart  and  West  1980, 
Shugart,  Mortlock,  Hopkins,  and  Burgess  1980,  Shugart  and 
Noble  1981,  Doyle  1982,  Doyle,  Shugart,  and  West  1982, 
Shugart  1984,  and  Pickett  and  White  1985)  have  been  utilized 
to  study  forested  ecosystems  around  the  world.   These  models 
have  various  gap  sizes  ranging  from  100m*2  to  833m"2. 


Spatial  Systems  and  Models 

A  spatial  predator-prey  insect  microcosm  was  used  by 
Huffaker  (1958)  to  study  two  species  of  mites.   The  prey 
mite  fed  on  oranges  while  the  predator  mite  fed  on  the  prey. 
In  one  set  of  experiments,  the  oranges  were  distributed  in  a 
10x12  grid  with  partial  barriers  between  the  oranges  and  one 
prey  placed  on  each  of  the  120  oranges.   Five  days  later  27 
predators  were  dispersed  on  the  oranges.  The  resulting  dy- 
namics in  populations  both  over  time  (8  months)  and  space 
are  shown  in  Figure  7.   In  other  experiments  with  oranges  in 
different  arrangements,  the  oscillatory  behavior  was  not 
seen.   Huffaker  concluded  that  the  predator-prey  oscillation 
would  only  occur  when  there  was  migration  from  the  outside 
or  a  sufficiently  complex  spatial  arrangement  of  prey  and 
barriers  to  allow  localized  growth  of  the  prey  followed  by 
consumption  by  the  predator. 
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In  high  altitude  balsam  fir  forests  in  the  northeastern 
United  States,  waves  of  tree  loss  and  regeneration  are 
thought  to  be  formed  by  an  interaction  of  the  prevailing 
wind  with  the  larger  mature  trees  that  are  exposed  along  the 
gap-wave  (Sprugel  and  Bormann  1981,  and  Sprugel  1984).   The 
wind  in  this  case  acts  to  organize  the  disturbance  cycle 
that  occurs  normally  in  this  type  of  forest  into  a  spatial 
wave  pattern  instead  of  randomly  occurring  patches. 

The  'ohi'a  dieback  phenomenon  in  the  rain  forests  of 
Hawaii  (Mueller-Dombois  1980)  is  a  case  of  localized  loss  of 
trees  in  the  forest  not  due  to  disease  or  insect  pest.   It 
was  postulated  that  the  effects  were  due  to  local  soil 
moisture  loss  arising  from  some  climate  instability.   Repro- 
duction of  the  'ohi'a  was  adequate  enough  to  regenerate  the 
forest  after  the  dieback,  thus  providing  a  way  for  this 
shade  intolerant  species  to  become  the  primary  canopy 
species  without  further  succession.   Climatic  variability 
was  thus  used  to  an  adaptive  advantage. 

Spatial  modelling  of  ecosystems  can  be  done  in  several 
different  ways.   By  using  a  model  based  on  the  FORET  simula- 
tion model  (Shugart  and  West  1977}  and  spatially  distribu- 
ting the  output  of  the  model  according  to  flooding  condi- 
tions and  hydroperiod,  Pearlstine,  McKellar  and  Kitchens 
(1985)  suggested  possible  species  changes  due  to  changes  in 
the  hydroperiod  caused  by  a  river  diversion  in  South 
Carolina.   In  this  case  the  number  of  individual  subcell 
models  was  kept  small  and  the  spatial  distribution  was  based 
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on  a  combination  of  terrain  relief,  hydrology,  and  cor- 
related output  from  the  simulation  model. 

Another  approach  to  spatial  modelling  is  to  divide  the 
area  into  individual  cells  with  a  representative  model  in 
each  cell  with  some  interaction  terms  among  the  individual 
cells.   This  is  the  approach  Costanza  (1979)  used  in  model- 
ling the  economic  development  of  South  Florida. 

Simulations  with  individual  models  for  each  cell  have 
certain  advantages,  because  the  interaction  of  neighboring 
cells  influences  the  outcome.  A   serious  disadvantage  where 
the  number  of  cells  is  large  is  the  immense  amount  of 
computer  time  required  for  the  simulations.   By  making  the 
cell  size  larger  this  can  be  avoided,  but  loss  of  spatial 
detail  occurs  as  the  cell  size  increases.   The  sub-cell 
distribution  modelling  technique  used  by  Pearlstine  et  al . 
(1985)  has  just  the  opposite  advantages  and  disadvantages. 
The  time  requirements  for  simulation  do  not  necessarily 
increase  as  the  area  of  cells  is  increased,  but  individual 
intercell  interactions  are  lost. 

Plan  of  Study 
Objectives 

This  study  of  energy  use  and  pattern  formation  with 
production  consumption  models  has  several  parts: 

First,  the  energetics  of  different  pathway  config- 
urations were  tested  using  a  series  of  miniraodels.   These 
models  were  manipulated  to  determine  the  energy  use  of 
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systems  with  different  production  and  consumption  kinetics 
and  different  combinations  of  components. 

Second,  a  generalize  production-consumption  minimodel 
calibrated  with  tropical  rainforest  data  was  used  to  study 
the  energetics  of  pulsing  behavior. 

Third,  spatial  pattern  formation  was  investigated  using 
the  pulsing  production-consumption  model  as  subunits  in  a 
spatially  distributed  format.   The  spatial  effects  and 
energy  implications  of  various  patterns  of  energy  inputs, 
edges,  and  lateral  connectivity  were  determined. 

These  spatial  simulations  included  several  types  of 
inter-block  exchange.   Hierarchical  relationships  are  rep- 
resented in  these  models  when  each  consumer  component  inter- 
acts with  more  than  one  producer  unit.   The  distribution  of 
gaps  developed  by  simulations  was  compared  with  gaps  in  the 
tropical  rainforest  in  Puerto  Rico. 

Finally,  insights  and  hypotheses  were  developed  about 
behavior  of  ecological  systems. 

Data  site:   Luquillo  Rainforest,  Puerto  Rico 

Data  from  the  Lower  Montane  Rainforest  in  the  Luquillo 
Mountains  of  Puerto  Rico  were  used  to  compare  some  of  the 
spatial  simulations  of  pulsing  and  patches.   Extensive 
studies  on  this  forest  were  published  previously  (Odum  and 
Pigeon,  1970) . 

Changes  in  structure  and  composition  of  a  plot  of 
tropical  rain  forest  near  El  Verde  in  Puerto  Rico  over  a 
period  of  30  years  were  reported  by  Crow  (1980).   Data 
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included  size  class  distributions  taken  in  1943,  1946,  1951 
and  1976  (Figure  8)-.   It  can  be  seen  that  there  is  a  shift 
over  time  in  the  different  size  classes.   The  peak  year  for 
the  0-8  cm  class  is  1946  while  the  peak  in  the  8-12  cm  class 
occurs  in  1951  and  the  peak  in  the  next  three  size  classes 
occurs  in  1976.   The  smallest  number  in  the  lower  two  clas- 
ses also  occurs  in  1976.   The  last  severe  hurricane  struck 
Puerto  Rico  in  1932,  and  this  movement  through  the  size 
classes  appears  to  be  the  growth  and  development  of  an  age 
class  of  trees  that  grew  back  after  the  hurricane.   The 
hurricane  in  this  case  acts  as  an  organizing  disturbance  to 
reset  succession  of  patches  on  a  large  scale. 

The  models  simulated  include  the  main  integrative  mech- 
anisms observed  in  ecosystems  for  coupling  production  and 
consumption  of  spatially  distributed  units.   Energy  use  of 
these  configurations  was  obtained  from  the  simulations  to 
test  the  hypothesis  that  commonly  observed  organizational 
designs  with  a  successional  regime  that  alternates  pro- 
duction and  consumption,  tend  to  maximize  system  power  in 
the  long  run. 


Figure  8.   Size  class  distribation  over  time  of  plot  of 
trees  in  tropical  forest  at  El  Verde  (Crow  1980) . 
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Distribution    of   trees   over  time 
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CHAPTER  2 
METHODS  AND  MODELS 

Ecosystem  concepts,  configurations,  and  models  were 
represented  with  energy  circuit  language  from  which  simula- 
tion programs  were  derived.  The  energy  circuit  language, 
developed  by  H.  T.  Odura  (Odum  1971,  Odum  and  Odura  1981  and 
Odum,  1983),  is  a  symbolic  language  for  modelling  ecosystems 
and  their  components.   Elements  of  storages,  flows,  and 
interactions  in  this  symbolic  language  keep  track  of  the 
laws  of  energy  conservation.   The  energy  diagrams  also  show 
the  correct  kinetic  interaction  between  parts  of  the  system. 
The  level  of  aggregation  or  disaggregation  that  is  needed  to 
understand  and  model  a  system  for  a  particular  purpose  can 
be  achieved  by  drawing  and  revising  diagrams  using  this 
energy  circuit  language.   A  diagram  of  most  of  the  important 
symbols  with  a  brief  description  of  each  is  presented  in 
Figure  9. 

One  of  the  benefits  of  using  the  energy  circuit 
language  is  that  it  is  possible  to  go  from  a  conceptual 
model  to  the  development  of  the  differential  equations 
needed  to  simulate  the  model  in  a  few  steps.   Each  of  the 
pathways  on  the  diagram  represents  a  flow  that  in  turn  can 
be  represented  by  terras  in  the  differential  equations  that 
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Figure  9.   Energy  circuit  langaage  symbols  (Odum  1983) 
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Energy  ctrcnit  A  pathway  whose  flow  is  propor- 

tional to  the  quantity  in  the  storage  or  somre  upstreazn. 

Source  Outsidt  source  of  ener^  delivehng  forces 

acconiing  to  a  progrHm  cootralled  ^m  out&ide;  a  forcing  func- 


TanJt  A  computment  of  energy  storage  within  the 

system  storing  a  quantity  bb  the  balance  of  InEows  and  outSnws^ 
a  Stat*  variables 

Heat  sink  Dispersion  of  potential  energy  into  heac 

that  accompanies  all  real  transformation  processes  and  storages; 
Iocs  of  pocantial  energy  from  further  use  by  the  system. 

InUraetion  Interactive  intersection  of  two  path- 

ways coiipied  to  producft  an  outflow  in  proportion  to  a  funcaon  of 
both;  control  actioaof  oneflowon  anothen  limiting  factor  action; 
work.  gate. 


Conaumer  Unit  that  transforms  energy 

quality,  stores- it.  and  feeds  it  bade  autocatalytically  to  improve 
inflow. 


Switching  action 
more  switching  actions. 


A  symbol  that  indicates  one  or 


r 


-    ) 

P- 

r 

Producer  Unit  that  collects  and  transforms 

low-quaiity- energy  under  control  interactions  of  high-quality 
flows. 


Self- limiting  energy  receiver  (Chapter  10).  A  unit  that  has  a 
self-lLraiting  output  when  mput  dnvea  are  high  because  there  is  a 
limiting  constant  quantity  of  material  reacting  on  a  circular 
pathway  within. 


Box  Miscellaneous  symbol  to  use  for  whatever  unit 

or  function  is  labeled- 


Constant-gain  amplifier  A  uni:  Lhat  delivers 

an  output  in  proportion  to  the  input  /  but  changed  by  a  constant 
factor  as  long  as  the  energy  source  S  is  suincient. 


Transection  A  unit  that  indicates  a  sale  of  goods 

or   services    (solid    line)    in    exchange    for    payment    of  money 
(dashed).  Price  is  shown  as  an  external  source. 
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describe  the  changes  in  storage  compartment  (tank)  values 
over  time. 

Simulation  Procedures  and  Programs 

The  majority  of  the  simulations  in  this  dissertation 
were  done  in  F0RTRAN-4-PLUS  on  a  Digital  Equipment  Corpora- 
tion (DEC)  POP  11/34  with  RSX-llM  operating  system.  The 
graphical  outputs  of  the  simulations  were  displayed  on  a  DEC 
VK-100  graphics  terminal  (General  Image  Generator  and  Inter- 
preter or  GIGI)  connected  to  a  Barco  color  monitor  and  DEC 
.LA-34  Decwriter.   The  GIGI  terminal  has  a  760x240  pixel 
resolution  and  can  display  up  to  eight  colors  on  a  color 
monitor.   In  order  to  facilitate  the  graphics  programming 
needed  in  my  simulation  models,  I  developed  a  set  of  FORTRAN 
subroutines  with  a  more  natural  calling  sequence  to  execute 
the  ReGIS  (Remote  Graphics  Instruction  Set  use  by  the  GIGI 
terminal)  commands  from  the  programs.  This  library  of  rou- 
tines (GGLIB)  is  listed  and  documented  in  the  Appendix. 

Some  of  the  goals  of  this  dissertation  were  to  examine 
the  structure  and  function  of  systems  in  time  and  space  and 
to  determine  how  variation  in  coefficients  may  affect  energy 
flows  and  storages  of  the  systems.   Graphical  display  pro- 
grams were  developed  to  project  a  simulated  3-D  surface  of 
the  output  of  various  state  variables  over  time  and  over  a 
range  of  input  conditions.   A  special  3-D  graphics  display 
program  was  written  to  display  the  output  of  these  model 
simulations  (program  PLOTZ,  Appendix), 
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The  spatial  models  are  broken  down  into  cells  that  show 
the  concentration  of  a  given  parameter  in  the  individual 
cell  as  a  color  block.  For  display  on  the  color  monitor  this 
provides  dramatic  views  of  the  model  changes  over  time  and 
space.   In  order  to  make  hardcopy  printouts  a  display  char- 
acter set  was  designed  so  the  density  of  the  dots  in  an 
individual  cell  was  correlated  to  the  color  of  the  cell. 
This  provided  a  way  of  screen-dumping  the  images  to  paper 
and  achieving  patterns  on  paper  that  were  similar  to  the 
ones  on  the  video  screen  (See  Appendix  for  a  listing  of  the 
character  set)  . 

Simulation  Models 
Minimodel  Tests 


First  a  group  of  minimodels  were  simulated  to  relate 
energy  use  to  basic  pathway  designs.  Then  spatial  models 
with  these  configurations  were  studied  for  energy  use  and 
pattern  formation. 

Three  path  minimodel 

In  order  to  understand  how  a  system  processes  variable 
energy  inputs,  builds  structure,  and  regulates  or  maximizes 
energy  flows,  a  simple  single  tank  model  was  simulated.   The 
model  is  similar  to  the  one  described  by  Odum  (1982)  that 
has  parallel  pathways  of  different  types  competing  for 
available  energy  (Figure  10).   The  model  has  a  flow  limited 
source  connected  to  a  single  storage  (tank)  by  three  dif- 
ferent pathways;  a  linear  pathway  (Jl) ,  an  autocatalytic 


Figure  10.   Three  pathway  model  jsed  to  test  effects  of 
various  energy  inputs  on  kinetic  mechanisms. 

Linear  input:  J1=K1*R 

Autocatalytio  input:      J2=K2*Q*R 
Quadratic  input:  J3=K3*Q*Q*R 

dQ=Jl+J2+J3-K4*Q 
R=J0-J1-K0*R*Q-K5*R*Q*Q 
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pathway  {J2)  ,  and  a  quadratic  pathway  (J3).   The  tank  has  a 
linear  drain. 

The  model  represents  a  system  that  can  change  its  use 
of  three  functional  pathways  to  get  energy.   The  linear 
pathway  represents  the  energy  flow  that  a  system   can  re- 
ceive without  any  feedback  in  this  pathway,  only  pathway 
resistance  to  the  flow.   Because  it  is  a  donor  controlled 
pathway,  the  system  has  no  control  on  the  flow.   Diffusion 
pathways  are  an  example  of  this  type  of  energy  flow.  The 
linear  pathway  is  very  efficient  because  it  takes  almost 
nothing  to  receive  the  energy. 

The  autocatalytic  pathway  has  a  feedback  from  the  sys- 
tem storage  for  interacting  with  an  energy  source  to  facili- 
tate the  capture  of  more  energy.  If  energy  is  available  to 
support  the  storage  this  pathway  may  lead  to  a  competitive 
advantage  over  the  linear  pathway.  The  efficiency  of  the 
autocatalytic  pathway  depends  on  the  energy  source,  the 
storage  and  the  pathway  coefficient.   A  pathway  of  this  type 
has  the  capability  of  capturing  more  available  energy. 

The  quadratic  pathway  has  a  self-stimulating  feedback 
(see  equation  on  Figure  10)  from  the  storage  to  capture 
available  energy.   Examples  of  cooperative  feeding  that  may 
fit  this  model  are  common  in  ecosystems  such  as  pack  hunting 
by  some  carnivores,  cell  and  organ  system  interactions  and 
the  cooperative  work  by  humans  in  developed  nations. 
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This  model  was  simulated  in  BASIC  (program  THREEPATH  in 
Appendix)  on  a  Digital  Equipment  Corporation  (DEC)  PDP  11/34 
using  a  DEC  VK-100  graptiics  terminal  (GIGI).  Measurements 
were  made  of  the  percent  of  the  input  power  used  while 
applying  various  levels  of  input  power  and  varying  the 
frequency  of  input  power.   Simulation  runs  were  also  made 
with  one  or  more  of  the  three  pathways  set  to  zero  to  deter- 
mine the  impact  of  the  various  pathways  on  the  overall 
system  behavior  and  power  utilization. 

In  conjunction  with  the  three  path  model  in  Figure  10, 
a  similar  model  with  the  same  inputs  but  with  additional 
higher  order  drain  pathways  was  simulated  to  determine  the 
effects  on  total  power  usage  (Figure  11).   In  any  system 
that  has  crowding  effects  or  high  storage  costs,  these 
drain  pathways  may  determine  how  the  system  processes  energy. 
The  model  has  a  linear  drain,  an  autocatalytic  drain  and  a 
quadratic  drain. 

The  basic  three  path  model  was  tested  for  the  effects 
of  size  and  turnover  time  on  the  percent  power  used  for 
various  power  inputs  by  varying  the  drain  coefficient  (K4  on 
Figure  10)  in  multiple  simulation  runs. 

The  percent  power  used  when  the  three  path  model  com- 
petes with  individual  storages  with  single  pathways  (Figure 
12)  was  also  simulated  to  see  how  the  various  pathways  may 
help  or  hinder  a  system.   The  competitors  are  individual 
tanks  with  single  pathways  corresponding  to  the  three  path- 
ways in  the  three  path  model. 


Figure  11.   Three  pathway  model  with  multiple  drain  path- 
ways. Used  to  test  effects  of  higher  order  drain  pathways  on 
threepath  model. 

Linear  input:  J1=K1*JR 

Autocatalytic  input:  J2=K2*Q*JR 

Quadratic  input:  J3=K3*Q*Q* JR 

Linear  drain:  J4=k4*Q 

Autocatalytic  drain:  J5=K5*Q*Q 

Quadratic  drain:  J5=K5*Q*Q*Q 

dQ=Jl+J2+J3-J4-J5-J6 
JR=J0-J1-K2'*J2-K3'*J3 
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Figure  12.   Three  pathway  model  with  individual  competing 
units  having  single  input  pathways  similar  to  combined 
model.   Coefficients  in  Appendix. 

Combination  tank: 

Linear  input:  J1=K1»JR 

Adtocatalytic  input:      J2=K2*Q*JR 
Quadratic  input:  .J3  =  K3*Q*Q* JR 

dQ=JH-J2  +  J3-K4*Q 

Single  tanks: 

Linear  input:  J1X=K1'*JR 

dQl=JlX-K4*Ql 

AQtocatalytic  input:      J2X=K2 ' *Q2* JR 
dQ2=J2X-K4*Q2 

Quadratic  input:  J3X=K3 ' *Q3*Q3* JR 

dQ3=J3X-K4*Q3 

JR=J0-J1-K2' *J2-K3'  *J3-J1X-K2' *J2X-K3' *J3X 


49 


50 


For  any  system  to  survive  over  the  long  term,  it  must 
fit  into  a  regime  of  disturbances  or  catastrophic  events 
from  sources  outside  its  own  boundaries.   The  system  must  be 
tuned  to  the  frequencies  of  those  systems  that  influence  it 
in  order  to  maximize  power  and  survive.   The  three  path 
model  was  simulated  with  various  frequencies  of  power  input 
to  see  how  the  various  pathways  process  power  at  different 
frequencies  and  amplitudes. 

Parallel  production-consumption  minimodel 

A  model  with  producers  in  parallel  was  used  to  study 
the  effects  of  competition  among  producers  (Figure  13).  The 
model  had  three  producers,  all  having  the  same  structure, 
with  one  aggregate  consumer  that  was  consuming  all  three  and 
feeding  back  as  a  multiplier  on  the  production  function  of 
each.   It  is  a  basic  predator-prey  model  with  competition 
among  the  different  producers,  along  with  feedback  control 
and  energy  constraints  in  the  form  of  a  flow  limited  source. 
Instead  of  having  combinations  of  pathways  that  can  vary, 
this  model  had  combinations  of  producers  that  could  vary. 

The  producers  had  different  turnover  times  and  coef- 
ficients so  that  Ql,  Q2,  and  Q3  represented  climax,  raid- 
successional  (shrub)  and  early  successional  (weed)  species. 
The  coefficient  of  consumption  (the  percent  of  each  producer 
the  consumer  eats  per  unit  time)  for  each  producer  was 
different.   The  weed  species  had  a  higher  value  than  the 
shrub  species,  which  was  higher  than  the  climax  species 


Figure  13.   Parallel  production-consurapt ion  model. 

Individual  rate  equations 
Rl  =  K1*Q1*JR*Q4 
R2  =  K2*Q2*JR*Q4 
R3  =  K3*Q3*JR*Q4 
R4  =  D1*Q1 
Rj  =  D2*Q2 
R6  =  D3*Q3 
R7  =  K7*Q1*Q4 
R3  =  K8*Q2*Q4 
R9  =  K9*Q3*Q4 

RIO  =  Fl* {K1*Q1*JR*Q4  +  K2*Q2*JR*Q4  +  K3*Q3*JR*Q) 
Rll  =  K0*(K7*Q1*Q4  +  K3*Q2*Q4  +  K9*Q3*Q4) 
R12  =  D4*Q4 
JR   =  J0/(1  +  L1*Q1*Q4  +  L2*Q2*Q4  +  L3*Q3*Q4) 


Rate  equations  for  state  variables 
dQl  =  Rl  -  R4  -  R7 
dQ2  =  R2  -  R5  -  R8 
dQ3  =  R3  -  R5  -  R9 
dQ4  =  Rll  -  R12  -  RIO 
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{Odum  1969).   A  list  of  coefficients  is  given  in  Appendix 
Table  3. 

Several  variations  of  this  basic  model  were  written  in 
FORTRAN  and  BASIC  computer  languages  and  simulated  on  both  a 
PDP  11/34  and  on  a  Heathkit  H8.   The  source  listing  for  the 
standard  parallel  production-consumption  model  {SUC10)  is 
presented  in  the  Appendix. 
Pulse  Model 

A  general  pulsing  ecosystem  model  (Figure  14)  was  de- 
signed to  test  various  hypotheses  about  energy  flows  and 
pulsing,  hierarchical  organization,  and  spatial  development 
of  ecosystems.   Some  of  the  structure  of  the  model  was 
derived  after  the  tests  of  the  threepath  model  and  the 
parallel  production-consumption  model.   The  model  had  many 
characteristics  of  ecosystems  such  as: 

1.  Flow  limited  resources  (representing  solar  based 
energy  resources) . 

2.  Nutrient  storage  within  the  boundaries  of  the  model. 

3.  Units  of  production,  consumption  and  storage. 

4.  Feedback  of  consumers  on  production  through  nutrient 
recycle . 

5.  Consumption  at  low  maintenance  rates  and  at  high 
pulsing  rates. 

6.  Production  through  a  fast  turnover  storage  into  a 
long  turnover  biomass  storage. 

The  basic  pulsing  ecosystem  model  was  tested  for  dif- 
ferent flow  rates,  initial  storages  and  energy  inputs.   From 
this,  a  baseline  understanding  of  the  dynamic  behavior  of 
the  model  and  energy  processing  capabilities  (as  percent 
power  used)  was  developed. 

The  pulse  model  (Figure  14)  was  similar  to  the  one  in 
Richardson  and  Odum  (1981)  with  some  changes  in  coefficients 
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and  flows  to  calibrate  it  to  a  tropical  rain  forest  ecosys- 
tem.  The  original  model  was  run  on  an  Electronics 
Associates  Incorporated  model  2000  Analog/Hybrid  computer. 
The  models  presented  in  this  dissertation  were  simulated  on 
a  DEC  POP  11/34.   The  multiple  simulations  of  the  pulse 
model  were  generated  with  a  version  of  the  program  that 
would  run  25  simulations  while  varying  a  coefficient  or 
initial  condition  over  those  25  runs  and  generate  data  files 
that  were  then  displayed  with  the  FORTRAN  program  PLOTZ  (See 
Appendix) .   The  source  listing  of  the  FORTRAN  pulse  program 
is  in  the  Appendix. 

The  pulse  model  was  calibrated  with  tropical  forest 
ecosystem  values  for  carbon  flows  and  storages  (Jordan  and 
Drewry  1969,  Odum  and  Pigeon  1970,  and  Brown,  Lugo,  Silander 
and  Liegel  1983).   The  energy  diagram  of  the  model  is  given 
in  Figure  14  and  the  equations,  coefficients  and  initial 
conditions  of  the  state  variables  are  given  in  Appendix 
Table  4. 

Pulse  Model  With  Prey-Predator  Sectors 

An  additional  higher  trophic  level  consumer  was  added 
to  the  pulsing  consumer  model  (Figure  14)  in  order  to  test 
the  relationship  of  turnover  time  and  hierarchical  matching 
of  consumers  (Figure  15).   The  extra  consumer  added  to  the 
model  had  the  same  structure  as  the  lower  level  pulsing 
consumer  (Q3),  with  both  linear  and  quadratic  pathways. 
This  model  was  tested  by  varying  the  turnover  time  of  the 
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highest  level  consumer  (Q5)  and  measuring  the  percent  power 
used  and  the  level  of  the  other  storages  in  the  system. 

Spatial  Models 

The  models  previously  discussed  were  time  domain  models 
with  no  spatial  effects.   However,  because  ecosystems  de- 
velop through  time  and  space  and  spatial  variations  can  be 
at  least  as  important  as  variations  in  time,  spatial  models 
were  developed  and  simulated  to  test  hypotheses  concerning 
spatial  development  of  ecosystems  such  as  energy  processing 
and  pattern  formation  and  hierarchical  control  of  pattern 
formation . 

The  basic  spatial  model  was  a  collection  of  subunits, 
each  one  a  pulsing  consumer  model  (Figure  14).   These  sub- 
units  were  organized  in  a  spatial  format.   When  this  simple 
model  was  simulated  in  a  spatial  format,  size  effects,  edge 
effects  and  the  consumer  range  of  influence  can  become 
important.   Intercell  interactions  between  individual  produ- 
cers, consumers,  nutrients,  and  energy  sources  may  be  im- 
portant in  energy  utilization  and  pattern  formation. 

Effects  of  edges  in  the  spatial  model  were  of  interest 
in  pattern  formation  and  energy  use.   Special  boundary  con- 
ditions were  defined  for  the  model  cells  along  the  edge. 
These  boundary  cells  were  manipulated  in  the  simulation 
model  in  order  to  study  the  effects  of  edges  on  energy  use 
and  pattern  formation.   The  boundary  cells  were  also  manip- 
ulated to  minimize  the  effect  of  edges  in  certain  runs  of 
the  model. 


60 


Any  ecosystem  can  be  divided  into  edge  and  non-edge 
(center)  parts.   The  amount  of  edge  in  an  ecosystem  is 
a  function  of  the  size  and  number  of  the  individual  patches 
within  it.   For  a  given  area,  as  the  number  of  subunits 
increases  the  percent  of  the  subunits  on  the  edge  decreases 
(  see  Figure  16)  . 

A  10x10  matrix  was  used  in  the  spatial  simulations, 
giving  36%  of  the  total  in  edge  cells  and  64%  in  non-edge 
cells.   This  size  model  was  chosen  to  reduce  the  edge  and 
yet  be  small  enough  to  simulate  in  a  reasonable  time.   Com- 
puter runs  for  this  model  lasted  approximately  3  hours  on  a 
PDP  11/34.   A  model  with  a  center  to  edge  ratio  of  10:1 
would  need  approximately  20  times  as  many  cells.   In  order 
to  test  the  effects  of  edges  on  the  model,  a  single  layer  of 
cells  was  added  around  the  outside  edges  of  the  10x10 
matrix,  giving  it  a  12x12  total  area  (Figure  17).  The  outer 
layer  was  not  acted  as  a  buffer  to  approximate  conditions  of 
an  edgeless  system. 

Arrangements  of  cells 

In  simulating  a  spatial  model,  many  arrangements  of 
cells  can  be  used.  The  simplest  form  used  was  a  linear  array 
with  cells  arranged  in  a  linear  ring.   For  two  dimensional 
models  the  cell  geometry  chosen  was  a  square.   This  was  done 
for  several  reasons: 


Figure  16.   Number  of  edge  and  center  cells  as  a  function  of 
total  number  of  cells  in  a  given  square  area. 
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Figure  17.   Cell  geometries  considered  for  spatial  models. 

(a)  Square  matrix  with  each  cell  having  4  side 
and  4  corner  neighbors.  Active  10x10  matrix 
embedded  in  a  12x12  matrix.   This  one  was 
chosen  for  the  spatial  simulations. 

(b)  Hexagonal  matrix  with  each  cell  having  6  sid^ 
neighbors. 
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1.  It  simplified  programming  the  model  because  two 
dimensional  arrays  in  FORTRAN  are  set  up  in  rows  and 
columns, 

2.  It  simplified  writing  the  graphics  routines  to  dis- 
play the  cells  on  a  graphics  terminal. 

3.  It  reduced  the  edge  effects  of  the  model. 

Ring  model 

A  modified  version  of  the  two  dimensional  spatial  model 
was  used  to  simulate  a  one  dimensional  case.   The  standard 
spatial  pulse  model  was  connected  head  to  tail  in  a  ring  of 
36  cells. 

Two  dimensional  models 

The  simplest  spatial  implementation  was  the  basic  pulse 
model  repeated  over  the  10x10  matrix  with  no  interactions 
between  individual  cells.   This  model  (program  DSPl)  was 
then  simulated  with  three  different  energy  forcing  functions: 

1.  The  energy  source  was  hierarchically  distributed 
(highest  energy  input  at  the  center  of  the  matrix). 

2.  The  energy  source  was  evenly  distributed. 

3.  The  energy  source  was  randomly  distributed. 

Energy  inputs  were  scaled  so  the  mean  input  over  the 
whole  matrix  could  be  held  constant  for  all  energy  types. 
Overall  energy  input  could  be  varied  to  test  pattern  devel- 
opment and  energy  use  with  various  energy  levels. 

Two  different  initial  conditions  were  tested.   A  suc- 
cessional  sequence  was  simulated  with  the  initial  values  of 
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stored  production  (bioraass,  Q2  in  Figure  14)  set  to  a  low 
level.  A  steady  state  configuration  was  also  used  in  which 
Q2  was  set  to  a  value  just  below  the  pulse  threshold.   The 
nutrient  tank  (Q4)  in  each  case  was  balanced  to  contain  the 
remainder  of  carbon  available  in  each  cell. 

This  model  tested  different  conditions  and  inputs. 

1.  Diffusion  was  allowed  between  nutrient  tanks  (Q4)  of 
each  subunit.  The  base  model  (DSPl)  allowed  nutrients 
to  diffuse  between  cells  at  various  diffusion  rates. 
The  outer  layer  of  non-reacting  cells  (see  Figure  17) 
had  constant  values  for  Q4  to  allow  tests  of  total 
diffusion  into  and  out  of  the  cell  matrix  (diffusion 
along  the  edges)  . 

2.  Diffusion  was  allowed  between  consumer  tanks  (Q3)  of 
each  subunit  (program  DSP1Q3)  .   The  outer  non-reacting 
cell  layer  was  set  to  a  constant  value  or  was  allowed 
to  float  (program  DSPIQZ)  at  the  average  of  the  inner 
10x10  matrix  to  simulate  a  continuous  sheet. 

Simulations  were  run  in  which  the  consumer  had  a  larger 
area  or  territory  than  the  producer.  A  model  variation 
(program  DSPIC)  was  tested  in  which  all  of  the  consumer 
tanks  were  clumped  into  one  tank  that  aggregated  consumption 
over  the  10x10  matrix  simultaneously.  This  version  also  had 
three  different  input  energy  patterns  available,  and  allowed 
diffusion  between  nutrient  (Q4)  tanks. 
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The  final  variation  was  a  model  with  production  com- 
partmentalized as  before  in  individual  cells  but  with  free 
roaming  consumers,  not  constrained  by  cell  boundaries.   One 
consumer  was  allowed  to  consume  and  move  about  the  matrix 
according  to  a  set  of  constraints.   When  the  consumer  grew 
above  a  preset  size,  it  was  split  into  two  equal  halves  and 
each  half  was  allowed  to  consume,  move  and  split  again.   An 
upper  limit  of  100  was  placed  on  the  total  number  of  con- 
sumers that  could  be  generated  during  the  run  (the  total  in 
the  10x10  matrix  of  the  previous  model  versions).   This 
model  also  had  three  different  energy  inputs  and  diffusion 
of  nutrients  (Q4 ) . 

Format  for  Spatial  Display  Graphs 

Data  from  the  spatial  pulsing  model  were  displayed 
using  the  format  shown  in  Figure  18.   The  spatial  distribu- 
tions of  the  producers  and  consumers  were  shown  at  various 
times  during  the  run  (usually  50  years  apart).   The  values 
of  producers  and  consumers  in  individual  cells  were  repre- 
sented by  the  density  of  dots  in  the  cell.   The  producer 
density  increment  was  2000  g/m'~2  with  a  range  of  3-16,000 
g/m"2  while  the  consumer  was  represented  by  an  increment  of 
50  g/m"2  and  a  range  of  0-400  g/m"2. 

Measurement  of  Hierarchies  at  El  Verde  Site 

In  order  to  compare  hierarchical  relationships  that 
were  generated  in  the  model  with  those  occurring  in  the 
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tropical  rain  forest  at  El  Verde,  several  measurements  were 
made  from  data  sets  from  the  tropical  rain  forest  study  at 
El  Verde  (1963-1967)  in  the  Luquillo  Mountains  of  Puerto 
Rico  (Odum  and  Pigeon  1970). 

A  data  set  (2048  samples)  characterizing  the  forest  at 
the  radiation  site  was  generated  by  the  U.  S.  Array  Corps  of 
Engineers  (Rushing  1970).   At  the  radiation  site,  every 
plant  1.8  m.  or  taller  was  enumerated  within  a  radius  of  30 
m.  from  the  center  of  the  site.   Each  plant  was  recorded 
with  the  species  name,  height,  diameter,  crown  diameter, 
exact  location,  and  various  other  parameters. 

Black  and  white  negatives  of  aerial  views  of  the  radia- 
tion site  (taken  November  1963  before  the  radiation  treat- 
ment) were  printed  as  8x10  inch  photographs.   Individual 
gaps  characterized  by  the  presence  of  Cecropia  peltata  (an 
early  successional  species)  were  digitized  from  the  photo- 
graphs using  a  personal  computer,  Coraplot  digitizer  and 
digitizing  program  written  especially  for  this  purpose 
(Measures  in  Appendix). 


CHAPTER  3 
RESULTS 


Simulation  of  Three  Path  Model 

Individual  Pathway  Tests 

The  amount  of  energy  flowing  through  each  of  the  path- 
ways in  the  three  path  model  (Figure  10)   depends  on    the 
total  energy  input  to  the  model.   As  input  power  (JO)  was 
increased  (Figure  19)  steady  state  flows  for  each  of  the 
pathways  changed.   Each  pathway  predominates  at  certain 
times.   The  linear  path  had  the  largest  power  flow  when 
input  power  was  low,  while  the  quadratic  pathway  had  the 
highest  flow  at  higher  power  inputs. 

When  input  power  was  increased  through  time  (Figure 
20),  there  was  no  steady  state,  but,  like  Figure  19  when 
power  increased,  the  energy  flow  shifted  from  the  linear 
pathway  to  the  autocatalytic  and  finally  to  the  quadratic 
path.  The  fraction  of  energy  remaining  (jr/JO)  also  de- 
creased over  time.   As  input  power  increased,  a  greater 
fraction  of  the  input  power  was  utilized. 

The  model  was  run  with  different  pathway  combinations 
(Figure  21)  and  with  various  power  inputs.   Each  curve  on 
the  graph  represents  a  steady  state  value  for  various  com- 
binations of  pathways  present  in  the  simulation.   Power  used 
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Figure  19.   Steady  state  power  utilization  o£  units  in  the 
three  path  model  (Figure  10)  as  a  function  of  input  power 
(JO)  . 
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Figure  21.   Steady  state  energy  flows  on  various  pathways 
and  combinations  of  pathways  in  the  three  path  model  (Figure 
10)  as  a  function  of  input  power  (JO) . 

Linear  pathway:  J1=K1*R 

Autocatalytic  pathway:      J2=K2*Q*R 
Quadratic  pathway:         J3=K3*Q*Q*R 
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at  any  given  input  was  highest  with  all  three  pathways  pres- 
ent.  For  any  combination  of  pathways  that  contained  the 
quadratic  path  (J1+J2+J3  or  J2+J3  or  J1+J3),  power  used 
increased  with  power  input  to  reach  the  same  asymptote  {>95% 
power  used).   A  slightly  lower  level  was  reached  for  path- 
ways dominated  by  the  autocatalytic  pathway  (J2  or  J2+J1). 
This  asymptote  was  approximately  90%  power  used  with  in- 
creasing power  input.   With  only  the  linear  pathway  enabled, 
no  change  occurred  in  percent  power  used  with  increasing 
power  . 

A  unique  situation  occurred  when  the  quadratic  pathway 
(J3)  existed  alone.   A  low  initial  storage  (Q)  did  not  pro- 
vide enough  feedback  on  the  J3  pathway  to  allow  growth. 
Percent  power  used  was  never  significant.   The  simulation 
with  only  J2  and  J2+J3  showed  zero  percent  power  used  at  low 
input  levels,  then  rose  quickly  at  higher  input  power. 

The  size  of  the  storage  (Q)  was  varied  to  see  the 
effects  on  energy  usage  (Figure  22).   This  was  achieved  by 
varying  the  depreciation  coefficient  (K4)  in  multiple  run 
while  increasing  power  input  in  the  three  path  model.   At 
high  values  of  K4  (fast  turnover  times),  increases  in  per- 
cent power  used  at  steady  state  with  increasing  power  were 
small.   With  decreasing  values  of  K4  (slower  turnover 
times) ,  percent  power  used  increased  for  the  initial  and 
final  values  of  input  power. 

The  addition  of  multiple  drains  with  different  struc- 
tures (Figure  11)  did  not  have  as  great  an  effect  on  the 


Figure  22   Simulation  of  three  path  model  in  Figure  10. 
Percent  power  used  as  a  function  of  energy  input  and  size  of 
drain  coefficient  (K4  varied  from  .02  to  2.0)'. 
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Figure  23.   Simulation  of  three  path  model  with  multiple 
drain  pathways  in  Figure  11.  Percent  power  used  as  a 
function  of  energy  input  (JO) . 

Linear  drain:  04=k4*Q 

Autocatalytic  drain:      D5=K5*Q*Q 
Quadratic  drain:  06=K6*Q*Q*Q 
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model  as  multiple  inflow  pathways.   The  percent  power  used 
was  lowest  when  all  combinations  of  drain  pathways  were 
enabled  (Figure  23).   Percent  power  used  increased  with 
increasing  input  power.   The  highest  value  for  percent  powerJ 
used  was  achieved  when  only  the  original  linear  drain  was  ' 
present.   Any  combination  with  the  linear  drain  used  less 
power  at  low  power  inputs  than  the  nonlinear  pathways  alone 
or  in  combination.   The  higher  order  drains  enabled  the 
system  to  draw  more  power  at  low  levels  than  when  combined 
with  linear   pathways.   This  effect  was  opposite  from  that 
with  input  pathways  at  very  low  power  where  the  nonlinear 
pathways  did  not  function  well  (see  Figure  21). 

The  effects  of  adding  competition  pathways  to  the  model 
(Figure  12)  can  be  seen  in  Figure  24.   In  this  case,  each  of 
the  competing  pathways  (single  tanks  Ql,  Q2,  and  Q3  with 
individual  pathways)  were  left  on  throughout  the  simula- 
tions.  Here  again  the  various  pathways  were  disabled  and 
simulations  run  with  varying  power  inputs.   The  results  were 
similar  in  some  ways  to  those  in  Figure  21  where  at  high 
power  inputs  the  percent  power  used  approached  one  of  two 
asymptotes.   The  greatest  percentage  of  power  utilization 
occurred  when  all  pathways  were  enabled  and  the  lowest  power 
utilization  occurred  when  only  J2  or  J3  were  enabled.  The 
addition  of  the  extra  competing  storages  increased  the  per- 
cent power  used  in  each  of  the  pathway  combinations  compared 
to  Figure  21.   These  extra  pathways  were  always  there  to  use 


Figure  24.  Simulation  of  three  path  competition  model  with 
various  pathways  enabled  (Figure  12) .  Percent  power  used  as 
a  function  of  energy  input  (JO) . 

Linear  pathway:  J1=K1*R 

Autocatalytic  pathway:      J2=K2*Q*R 
Quadratic  pathway:  J3=K3*Q*Q*R 
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whatever   power    maybe    left  over    (particularly  the   linear 
path)  . 

Frequency  Studies 

The   basic    three   path   model    (Figure    10)    was   also   used    to 
test   the   effects   of   different   frequencies   of    input   power   on 
the    model   at    three    different  power    levels.        At   the   lowest 
power    level    (J0=500,    Figure    25)    the    differences   between 
pathways    in   percent   power   used   was    the   greatest.      The   great- 
est   frequency   response   occurred   at    low    frequencies.      The 
frequency   response    was    flat    with   only   the    linear   path   en- 
abled.     When   all  pathways   were  present,    the  percent  power 
used  was  highest  with  a  peak  at  approximately  2  cycles.     A 
peak  of  power   utilization   also  occurred   with    the   combin- 
ations of  J1+J2  and  J1+J3.     The  pathways  that  showed  a  min- 
imum  in   the   frequency   response   were   composed   of   J2+J3    (the 
two   nonlinear   pathways   combined)    and   J2.      The   quadratic 
pathway  alone   did  nothing   since   no  power    was   used    (compare 
with    Figure    21). 

When    the    input   power    was    increased    to   2000    (Figure    26), 
the    linear   pathway   showed  no  change   in   output    with   change    in 
frequency   and    the   quadratic   pathway   had   no   output.      The 
combination  of  J1+J2   here   again   had  a   slight   maximum   at 
about   2  cycles   while  J2  alone  had  a  maximum  at  zero  cycles. 
The    combination   of    all   of    the    pathways    (J1+J2+J3)    and    J1+J3 
had   a   slight   minimum   of  power   utilization   at  about   8   cycles, 
while   the  combination  of  J2+J3   showed  a   slight  minimum   at 
about    2    cycles. 


Figure  25.  Simulation  of  the  three  path  model  in  Figure 
10.  Percent  power  used  as  a  function  of  frequency  of  the 
input  power  (J0=500) . 


Linear  pathway: 
Autocatalytic  pathway: 
Quadratic  pathway: 


J1=K1*R 

J2=K2*Q*R 

J3=K3*Q*Q*R 
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Figure  26.   Simulation  of  the  three  path  model  in  Figure 
10.   Percent  power  used  as  a  function  of  frequency  of  the 
input  power  (J0=2000). 

Linear  pathway:  J1=K1*R 

Autocatalytic  pathway:      J2=K2*Q*R 
Quadratic  pathway:  J3=K3*Q*Q*R 
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When  the  input  power  was  raised  to  10000   (Figure  27) 
the  percent  power   used   went   up   for   all   combinations   of 
pathways   except    the    linear   path.      The   quadratic  pathway   was 
operational  at   this   high  power    level   but   with   a    significant 
minimum  at   2  cycles  per   run.      Other   combinations   had   small 
minima  and  maxima  that  are  hard  to  see  at  the  scale  of  this 
graph. 

The   response   of    the    model   to   various   frequencies   and 
power    input   is   shown   in   Table    1.      Simulation    runs    with 
pathway  J1+J2   had  a   maximum    in  percent  power   used   at  all 
three  power    inputs   while   the   combination  of   J2+J3   had   a 
minimum    in   percent   power   used  at  all   three   power    inputs. 
The   combination   of   all  pathways    (J1+J2+J3)    has   a   peak   of 
maximum  percent  power   utilization  at   low  power   and   low 
frequency    input.      At   higher   power    levels   percent   power    util- 
ization   (with   all  three  pathways  enabled)    was    lower   with 
some    shifting    in   the    frequency   at    which    this   occurs. 

Simulation  of  Parallel  Production-Consumption   Model 

Single  Run  Simulations 

The   parallel  production   model   showed  a   successional 
pattern   with   the    initial  dominant   species    (Q3,    with    the 
fastest    turnover)    growing    up,    then   declining    as   Q2    became 
the   dominant   species   and   finally   Ql    (with    the    slowest    turn- 
over)   reached   a   maximum   and    then   dropped  back   to  a   slightly 
lower    steady    state    (Figure    28).       The    consumer    (Q4,     with    the 


Figure  27.   Simulation  o£  the  three  path  model  in  Figure 
10.   Percent  power  used  as  a  function  of  frequency  of  the 
input  power  ( J0=10000) . 

Linear  pathway:  J1=K1*R 

Autocatalytic  pathway:      J2=K2*Q*R 
Quadratic  pathway:         J3=K3*Q*Q*R 
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Table  1.   Frequency  response  (rainimuras  and  maximums)  of 
three  path  model  (Figure  10)  with  varying  input  power. 


Pathway 
combination 


INPUT  POWER 
J0=500  J0=2000         J0=10000 


J1+J2+J3  (All)  MAX  (2)  MIN  (8)  MIN  (3) 

J2+J3  MIN  (2)  MIN  (2)  MIN  (2) 

J1+J3  MAX  (2)  MIN  (8)  MIN  (3) 

J1+J2  .  MAX  (2)  MAX  (2)  :^1AX  (2) 

Jl  N/R  N/R  N/R 

J2  MIN  (2)  MAX  (0)  MAX  (0) 

J3  N/0  N/0  MIN  (2) 


Numbers  in  parenthesis  are  the  frequencies  at  which  the 
maximum  or  minimum  occurs. 

N/R  signifies  there  was  no    frequency  response  for  this  set 
of  pathways. 

N/0   signifies  there  was  no  power  uses  at  these  inputs 
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longest  turnover  time  per  unit)  also  rose  to  a  steady  state 
value.   During  this  time,  the  productivity  climbed  to  a 
local  maxima,  then  dropped  slightly,  finally  climbing  to  a 
slightly  higher  steady  state.   The  percent  power  used  for 
the  whole  run  was  95.5%. 

In  order  to  test  the  role  of  each  of  the  producers 
early  in  the  simulation,  a  series  of  runs  were  made  with  the 
initial  condition  of  one  of  the  producer  species  set  to  zero 
(Figure  29  a,b,c).   With  no  initial  climax  species  (Ql) 
present  (Figure  29a)  the  shrub  species  (Q2)  became  dominant 
in  the  final  steady  state.   The  percent  power  used  for  the 
run  was  94.7%,  slightly  less  than  the  base  run  configura- 
tion.  This  configuration  did  not  support  as  high  a  level  of 
consumer  (Q4)  compared  to  the  base  model  run  (75.2  vs. 
90.8)  . 

When  the  shrub  species  (Q2)  was  absent  (Figure  29b), 
the  percent  power  used  for  the  run  and  steady  state  values 
for  the  consumers  were  similar  to  the  base  run.   Without  the 
shrub  species  present  to  compete  during  the  middle  period, 
the  final  climax  species  (Ql)  peal<ed  earlier  and  higher  than 
in  the  base  run. 

When  the  weed  species  (Q3)  was  initially  absent  (Figure 
29c),  the  system  was  not  self  sustaining.   The  primary 
reason  was  that  during  the  early  part  of  the  simulation,  the 
consumer  (Q4)  was  dependent  on  the  weed  species  (Q3).   With 
no  Q3  present,  the  consumer  crashed  very  quiclily.   The  whole 
system  then  crashed  because  the  consumer  feeds  back  in 


Figure  29   Simulation  of  the  parallel  production-consamption 
model  in  Figure  13.   See  Figure  28  for  legend  and  ordinate 
scale. 

(a)  Simulation  run  with  initial  value  of  climax 
species  (Ql)  set  equal  to  zero. 

(b)  Simulation  run  with  initial  value  of  interme- 
diate species  (Q2)  set  equal  to  zero. 

(c)  Simulation  run  with  initial  value  of  weed 
species  {Q3)  set  equal  to  zero. 
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the  production  function  of  all  of  the  producers  in  the 
system. 

When  the  model  was  simulated  with  no  initial  consumer 
(Q4) ,  it  crashed  even  faster  (not  shown)  than  in  Figure  29c 
because  of  the  feedbacks  in  the  model  from  the  consumer  to 
the  producers. 

Multiple  Run  Simulations 

The  behavior  of  the  parallel  production  model  with 
varying  input  power  is  seen  in  Figure  30a-f.   In  this  set  of 
runs  the  base  model  was  run  for  100  time  units.  For  each 
successive  run,  the  input  power  (JO)  was  increased,  varying 
from  50  to  300.   As  the  energy  input  increased,  the  peaks  of 
the  producers  were  higher  (Q1-Q3) ,  with  Q3  (the  weed 
species)  showing  the  most  change  in  amplitude  (Figure  30c). 
The  climax  species  (Ql)  peaked  sooner  as  the  input  power  in- 
creased . 

The  simulation  of  succession  to  a  climax  was  thus 
speeded  up  by  increasing  the  energy  input  at  lower  levels, 
but  at  higher  levels  the  increase  in  energy  had  little 
effect  on  the  transition  to  dominance  of  the  climax  species. 

The  effect  of  increasing  energy  input  was  also  seen  in 
the  level  of  the  consumer  (Q4  in  Figure  30d)  .   With  in- 
creasing power,  the  consumer  was  maintained  at  a  propor- 
tionately higher  steady  state. 

For  this  set  of  simulations,  as  the  input  power  in- 
creased, the  percent  power  used  increased  asymptotically 
(Figure  31).  There  was  a  diminishing  return  on  the  input 


Figure  30.   Simulation  of  the  parallel  production- 
consumption  model  in  Figure  13.   Multiple  simulations  o£  the 
model  with  available  power  increasing  Erora  50  to  300. 

(a)  Climax  species  (Ql) 

(b)  Intermediate  producer  species  (Q2) 

(c)  Weed  species  {Q3) 

(d)  Consumer  species  (Q4) 

(e)  Percent  power  used  (JO-Jr)/JO 

(f)  Total  biomass 
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Figure  31.   Simulation  of  the  parallel  production- 
consumption  model  in  Figure  13.   Multiple  simulations   of 
t'ne  model  with  percent  power  used  for  entire  run  vs  input 
power.   See  Figure  RS3a-f. 
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power  as  the  effect  was  greater  at  low  power  than  it  was  at 
higher  levels  of  power. 

When  the  input  power  was  varied  as  in  the  previous 
example  (50  to  300)  but  the  initial  condition  of  the  con- 
sumer was  started  at  a  higher  level  (Q4INIT=50,  lOx  base  run 
value)  the  results  were  similar  to  the  previous  run  but 
damped  (Figure  32a-f).   The  shift  in  time  of  the  peak  of  the 
climax  species  (Ql)  was  less  than  before  and  the  amplitudes 
of  the  initial  peaks  of  Q2  and  Q3  were  less.   Percent  power 
used  per  time  increment  also  was  higher  in  the  earlier 
stages  of  this  run  compared  to  the  previous  run  (compare 
Figure  32e  with  30e).   With  higher  initial  levels  of  the 
consumer,  the  model  generated  more  power  earlier  through  the 
feedback  of  the  consumer  on  the  producers. 

When  the  input  power  was  held  constant  (J0=100)  and  the 
initial  condition  of  the  consumer  (Q4)  varied,  the  model 
displayed  two  different  behaviors  (Figure  33a-f).   With  few 
consumers  initially,  the  system  crashed,  unable  to  proceed 
through  the  normal  growth  sequence.   When  the  initial  quant- 
ity of  the  consumers  (Q4)  was  above  a  critical  level,  the 
system  grew  and  went  through  a  normal  growth  sequence.   A 
sharp  transition  occurred  in  the  percent  power  used  as 
Q4INIT  was  increased  (Figure  34). 

Because  the  consumer  (Q4)  was  feeding  back  as  a  multi- 
plier to  the  producers,  some  minimum  critical  value  must 
exist  for  the  consumer  population  to  stabilize  this  model. 


Figure  32.   Siraulatioa  of  the  parallel  prodaction- 
consaraption  model  in  Figure  13.   Run  with  available  power 
increasing  from  50  to  300  and  the  initial  value  of  the 
consumer  (Q4)  equal  to  50  ( lOx  base  run  in  Figure  28). 

(a)  Climax  species  (Ql) 

(b)  Intermediate  producer  species  (Q2) 

(c)  Weed  species  (Q3) 

(d)  Consumer  species  (Q4) 

(e)  Percent  power  used  {JO-Jr)/JO 

(f)  Total  bioraass 
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Figure  33.   Simulation  of  the  parallel  production- 
consumption  model  in  Figure  13.   Multiple  simulations  of  the 
model  with  available  power  held  constant  {J0=100,  base  run 
value)  and  the  initial  value  of  the  consumer  (Q4)  varied 
from  1  to  6. 

(a)  Climax  species  (Ql) 

(b)  Intermediate  producer  species  (Q2) 

(c)  Weed  species  (Q3) 

(d)  Consumer  species  (Q4) 

(e)  Percent  power  used  (JO-Jr)/JO 

(f)  Total  bioraass 
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Figure  34.   Simulation  of  the  parallel  production- 
consumption  model  in  Figure  13.   Total  percent  power  used 
for  entire  run  as  a  function  of  the  initial  value  of  the 
consumer  (Q4) .   This  represents  a  cross  section  of  Figure 
33e. 
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Either  immigration  or  a  temporary  auxiliary  support  system 
is  necessary  to  start  a  system  of  this  class. 

Similarly,  when  the  input  power  was  held  constant 
(J0=100)  and  the  initial  condition  of  the  weed  species  (Q3) 
was  varied  (Figure  35a-f)  ,  the  system  crashed  at  low  levels 
of  Q3,  but  at  higher  levels  it  was  stable  (see  Figure  29c 
for  a  single  run  with  Q3=0). 

The  system  response  was  different  with  changes  in  the 
initial  conditions  of  Ql  and  Q2  (refer  to  Figures  29a  and 
29b)  because  the  consumer  was  not  as  dependent  upon  them  for 
its  survival  early  in  the  simulation. 

Initial  Conditions  and  Total  Energy  Use 

The  behavior  of  the  parallel  production-consumption 
model  with  different  initial  conditions  for  the  state  var- 
iables (Ql,  Q2,  Q3,  and  Q4)  and  input  power  was  tested.   In 
this  set  of  simulations,  the  total  percent  power  used  was 
measured  for  each  simulation  run  while  varying  the  input 
power  and  the  initial  condition  of  the  state  variables  one 
at  a  time  (Figure  36). 

In  all  four  cases  when  JO  was  low,  the  model  was  unable 
to  utilize  the  energy  available  to  it.  When  the  input  power 
was  above  a  certain  point  then  the  model  was  able  to  utilize 
the  input  energy  with  two  exceptions.  When  Q3  (weed 
species)  was  very  low,  the  percent  power  used  rose  to  a 
plateau  then  fell  when  the  input  energy  went  above  a  certain 
level.   The  model  was  unstable  under  these  conditions. 


Figure  35.   Simulation  of  the  parallel  production- 
consumption  model  in  Figure  13.  The  initial  value  of  weed 
species  (Q3)  was  varied  from  0  to  .5  and  input  power  was 
held  constant  (J0=100,  base  run  value). 

(a)  Climax  species  (Ql) 

(b)  Intermediate  producer  species  (Q2) 

(c)  Weed  species  (Q3) 

(d)  Consumer  species  (Q4) 

(e)  Percent  power  used  (JO-Jr)/JO 

(f)  Total  biomass 
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Figure  36.  Steady  state  values  of  percent  power  used  as  a 
function  oE  input  energy  and  state  variable  initial 
conditions  for  multiple  simulation  runs  of  parallel 
production-consumption  model  (Figure  13) . 

(a)  Vary  input  energy  and  Ql  (Climax  species) 

(b)  Vary  input  energy  and  Q2  (Intermediate  producer) 

(c)  Vary  input  energy  and  Q3  (Weed  Species) 

(d)  Vary  input  energy  and  Q4  (Consumer) 


115 


100 


300 


b. 


300 


JO 


117 

When  Q4  was  below  a  certain  threshold  the  system  could 
not  be  sustained  regardless  of  the  input  energy.   After  an 
initial  threshold  level  of  consumers  was  reached,  the  system 
was  stable,  similar  to  that  described  above  for  Q3.   Since 
Q4  has  a  direct  feedback  on  Ql,  Q2,  and  Q3,  the  interaction 
of  these  in  the  production  term  can  determine  whether  or  not 
the  system  was  stable.   If  the  value  of  Q4  was  too  low  then 
there  was  little  production  and  the  system  crashed. 

Simulation  of  the  Pulse  model 

Single  Run  Simulations 

A  simulation  of  the  base  run  pulse  model  (Figure  14)  is 
shown  in  Figure  37a.   As  Q2  increased,  the  available  carbon 
or  nutrient  carbon  tank  {Q4)  decreased  proportionately.   As 
the  stored  biomass  increased  there  was  a  threshold  level  at 
which  the  consumer  (Q3)  began  to  grow  rapidly  and  pulsed. 
This  pulse  consumed  Q2  and  released  the  carbon  back  into  the 
available  carbon  pool  (Q4).   The  threshold  of  pulsing  was 
dependent  on  the  level  of  both  Q2  and  Q3.   The  level  of  Q3 
before  the  pulse  was,  however,  directly  related  to  the  level 
of  Q2  and  the  input  diffusion  pathway.   After  the  pulse,  the 
consumer  (Q3)  decayed  back  to  a  low  level. 

The  cycle  repeats  itself  at  a  frequency  of  approximate- 
ly 325  years.   The  power  used  varied  during  the  simulation 
with  the  highest  rate  occurring  shortly  after  the  pulse, 
when  the  nutrients  have  been  concentrated  in  Q4  as  available 
carbon . 


Figure  37.   Simulation  for  pulse  model  (Figure  14)  with 
base  run  coefficients  (See  Appendix). 

(a)  Base  run  of  model. 

(b)  Input  energy  one-half  of  base  run. 

(c)  Input  energy  two  times  the  base  run. 

Legend: 

PPQ  =  Percent  power  used 

Ql  =  Production  unit 

Q2   =  Stored  biomass 

Q3   =  Pulse  consumer 

Q4   =  Nutrient  storage 
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If  the  input  energy  was  less  (J0=50,  half  of  the  base 
run)  then  the  pulse  came  at  a  later  time  (Figure  37b)  and 
the  frequency  of  pulsing  had  a  longer  period.   The  production 
was  lower  and  the  stored  biomass  {Q2)  took  longer  to  reach 
the  level  that  would  trigger  the  pulse  in  the  consumer  (Q3). 

When  the  input  energy  was  raised  to  twice  the  level  of 
the  base  run  (J0=200),  the  consumer  pulsed  only  one  time 
(Figure  37c)  and  then  remained  at  a  low  level  instead  of 
decaying  away  entirely  as  in  the  base  run.  With  a  low  level 
of  consumer,  the  stored  biomass  was  not  able  to  build  up  and 
remained  at  a  lower  steady-state  level. 

The  total  power  used  for  each  of  these  runs  was  related 
to  the  input  power.  As  the  input  power  went  up,  the  percent 
power  used  also  went  up  from  93.3  at  50%,  to  96.5  at  100% 
and  98.3  at  200%. 

The  quadratic  pathway  between  the  stored  biomass  (Q2) 
and  the  consumer  (Q3)  was  responsible  for  the  pulsing   much 
as  the  autocatalytic  pathway  of  a  Lotka-Volterra  model  is 
responsible  for  its  oscillating  limit-cycle  behavior.   With 
only  the  linear  path  between  Q2  and  Q3,  the  behavior  was  not 
pulsing  or  oscillatory  (Figure  38).   The  stored  biomass  grew 
while  the  nutrients  were  used  up.   In  this  time  frame  (760 
years),  the  values  did  not  reach  a  steady  state  and  93%  of 
the  available  power  was  used,   when  simulated  for  2000  time 
units  (Figure  38b) ,  the  percent  power  used  dropped  off  to  a 
low  steady  state  value.   The  system  became  nutrient  limited 


Figure  33.   Simulation  of  pulse  model  (Figure  14)  without 
a  quadratic  pathway  {K7,  K8,  K9  =  0.0). 

(a)  Simulation  for  760  years. 

(b)  Simulation  for  2000  years. 

Legend: 

PPU  =  Percent  power  used 

Ql  =  Production  unit 

Q2   =  Stored  bioraass 

Q3   =  Pulse  consumer 

Q4   =  Nutrient  storage 
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because  most  of  the  nutrients  were  tied  up  in  the  stored 
biomass. 

When  the  pulse  model  was  run  without  feedbacks  into  Q4 
(pathways  R6  and  R8  cut  off)  the  model  continued  to  pulse 
but  began  to  decline  (Figure  39).   The  percent  power  used 
dropped  as  the  level  of  Q4  dropped  until  one  final  pulse  and 
then  everything  decayed  to  a  low  steady  state  condition. 

Multiple-run  Simulations 

When  the  input  power  was  increased,  the  result  was  most 
noticeable  on  the  stored  biomass  (Q2)  and  the  consumer  (Q3, 
Figure  40).  At  low  values  of  JO  there  was  no  pulsing  within 
the  time  frame  of  the  simulation  (760  years).   As  JO  was 
increased,  the  pulsing  began  as  a  result  of  the  stored 
biomass  (Q2)  increasing  to  a  threshold  level  at  which  Q3 
pulsed  and  consumed  the  stored  biomass  (Q2).   As  JO  was 
further  increased,  the  pulsing  frequency  increased.   At  high 
levels  of  JO  the  first  pulse  decayed  and  the  system  switched 
to  a  steady  state  with  Q2  being  maintained  at  a  low  level 
(see  Figure  37c  for  example).   The  total  power  used  (Figure 
40e)  increased  linearly  as  JO  increased  with  small  fluc- 
tuations over  time  due  to  the  pulsing  of  Q3.   The  percent 
power  used  (Figure  40f)  was  less  than  80%  for  low  values  of 
JO  then  rose  rapidly  through  the  pulsing  and  leveled  off   as 
JO  approached  250.   The  percent  power  used  was  reduced  by 
the  initial  consumption  but  returned  to  a  maximum  after  the 
pulse.   The  percent  power  used  increased  as  the  available 


Figure   39.      Simulation   of  pulse   model    (Figure   14)    without 
feedbacks    into    Q4    (K6,    K8    =    0.0) 

(a)  Simulation  for  760  years. 

(b)  Simulation  for  2000  years. 

Legend; 

PPU  =  Percent  power  used 
Ql   =  Production  unit 
Q2   =  Stored  biomass 
Q3   =  Pulse  consumer 
Q4   =  Nutrient  storage 
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Figure  40.   Multi-run  simulation  of  the  pulse  model  (Figure 
14)  with  variation  in  input  energy.  {JO  varied  from  0  to 
250)  . 

(a)  Production  unit  (Ql) . 

(b)  Stored  biomass  (Q2). 

(c)  Pulse  consumer  (Q3) . 

(d)  Nutrient  storage  (Q4) . 

(e)  Power  used  (JO-Jr) 

(f)  Percent  power  used  100* (JO-Jr) / (JO) 
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power  was  increased  (similar  to  three  path  models  seen 
earlier)  with  local  maxima  immediately  after  the  pulse. 

The  total  amount  of  nutrients  in  the  system  also  had  an 
important  effect  on  the  behavior  of  the  model  (Figure  41a- 
f )  .   At  higher  initial  levels  of  Q4  there  was  little  effect 
on  the  model.   At  these  higher  ranges,  the  model  was  no 
longer  nutrient  limited  but  was  energy  limited.   At  low 
values  for  the  initial  concentration  of  Q4  the  pulsing 
greatly  affected  the  labile  production  (Ql),  the  stored 
biomass  (Q2)  and  the  pulsing  consumer  (Q3)  .   At  the  lowest 
level  of  Q4,  there  was  no  pulsing,  Q2  remained  at  a  low 
steady  state  value,  and  Q3  also  remained  at  a  low  steady 
state  value.   There  was  a  small  shift  in  the  pulsing  fre- 
quency at  the  lowest  initial  levels  of  Q4  but  no  frequency 
shift  at  the  higher  levels.   The  power  used  was  greatly 
affected  at  low  initial  levels  of  Q4  but  rose  only  slightly 
at  higher  values  of  Q4.   For  the  same  amount  of  change  in 
Q4,  the  variability  of  the  power  used  was  greater  when  Q4 
was  small  than  when  Q4  was  high.   However,  the  percentage 
change  was  greater  in  the  beginning  than  at  the  end. 

The  turnover  time  of  the  pulsing  consumer  affected  the 
behavior  of  the  system  and  use  of  power  (Figure  42a- f).  The 
pulse  model  was  simulated  with  the  value  of  the  drain  coef- 
ficient (K12)  of  the  consumer  (Q3)  varied  with  each  run.   As 
the  turnover  time  increased,  the  frequency  of  pulsing  shift- 
ed to  a  shorter  period  with  the  amplitude  decreasing  until 
there  is  no  pulse  at  all  but  a  continually  rising  consumer. 


Figure  41.  Multi-run  simulation  of  pulse  model  (Figure 
14)  with  variation  in  total  carbon  in  model.  (Q4  varied 
from  2000  gC/m2  to  100,000  gC/ra2. 

(a)  Production  unit  (Ql)  . 

(b)  Stored  bioraass  (Q2) . 

(c)  Pulse  consumer  (Q3) . 

(d)  Nutrient  storage  CQ4) . 

(e)  Power  used  (JO-Jr) 

(I)       Percent  power  used  100* (JO-Jr) / (JO) 
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Figure  42.   Multi-run  simulation  of  pulse  model  (Figure 
14)  with  variation  is  turnover  time  of  pulsing  consumer. 
(K12  varied  from  .01  to  .5). 

(a)  Production  unit  (Ql) . 

(b)  Stored  biomass  (Q2) . 

(c)  Pulse  consumer  {Q3) . 

(d)  Nutrient  storage  (Q4). 

(e)  Power  used  (JO-Jr) 

(f)  Percent  power  used  100* (JO-Jr )/ (JO) 
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This  implies  there  is  a  'window'  of  size  for  the  consumer  to 
pulse . 

Changing  the  rate  constant  (K9)  of  the  quadratic  path- 
way caused  the  pulsing  consumer  to  change  frequency,  in- 
creasing the  frequency  of  pulsing  with  an  increasing  coef- 
ficient value  (Figure  43).   There  was  a  point  in  this  set  of 
simulations  where  the  pulsing  ceases  but  in  this  case  the 
size  of  the  consumer  remains  small.   When  the  quadratic 
pathway  became  dominant  at  low  consumer  levels,  the  system 
did  not  pulse  but  completely  consumed  the  stored  bioraass 
storage  (Q2). 

When  simulated  without  the  quadratic  pathway  and  chang- 
ing the  coefficient  of  the  linear  pathway  (Kll) ,  the  model 
did  not  pulse,  the  consumer  (Q3)  remained  at  a  low  level  and 
the  stored  biomass  (Q2)  built  up  (Figure  44,  compare  to 
single  run  Figure  38).   As  the  linear  pathway  increased, 
there  was  a  slight  increase  in  the  consumer  (Q3)  with  less 
of  a  build-up  in  the  stored  bioraass  (Q2).   In  all  cases, 
through  time  the  power  use  and  percent  power  used  dropped 
off. 

Simulation  of  Pulse  Model  with  Prey-Predator  Sectors 

Simulation  of  the  pulse  model  with  an  additional  prey- 
predator  sector  (Figure  15)  investigated  how  turnover  time 
is  related  to  hierarchical  consumers  (Figure  45).   With  a 
drain  coefficient  on  Q5  the  same  as  or  larger  than  that  of 
the  normal  pulsing  consumer  (K17  =  0.05  or  0.5),  the  effect 


Figure  43.  Multi-run  simulation  of  pulse  model  (Figure 
14)  with  variation  in  quadratic  pathway  (K9  varied  from 
0.5E-6  to   0.53E-5  with  K7  and  K8  varied  proportionately) 

(a)  Production  unit  (Ql). 

(b)  Stored  bioraass  (Q2) . 

(c)  Pulse  consumer  (Q3). 

(d)  Nutrient  storage  (Q4). 
(a)  Power  used  (JO-Jr) 

(f)   Percent  power  used  100* (JO-Jr) / (JO) 
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Figure  44.   Multi-run  simulation  of  pulse  model  (Figure 
14)  with  variation  in  linear  pathway  (Kll  varied  from  0.0 
to  0.12E-2  and  K5  and  K5  varied  proportionately)  with  quad- 
ratic pathway  held  at  zero. 

(a)  Production  unit  (Ql) . 

(b)  Stored  bioraass  (Q2) . 

(c)  Pulse  consumer  {Q3) . 

(d)  Nutrient  storage  (Q4) . 

(e)  Power  used  (JO-Jr) 

(f)  Percent  power  used  100* (JO-Jr) / (JO) 
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Figure  45.   Simulation  of  pulse  model  with  prey-predator 
sectors  (Figure  15) . 

(a)  Simulation  with  turn-over  time  of  higher  level 
pulsing  consumer  (Q5)  set  equal  to  lower  level 
pulsing  consumer  (Q3) . 

(b)  Simulation  with  turn-over  time  of  higher  level 
pulsing  consumer  (Q5)  set  to  ten  times  longer  than 
the  turn-over  time  of  lower  level  oulsing  consumer 
(03)  . 

(c)  Simulation  with  turn-over  time  of  higher  level 
pulsing  consumer  (Q5)  set  to  one  hundred  times 
longer  than  the   turn-over  time  of  the  lower  level 
puling  consumer  {Q3) . 
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was  hardly  detectable  in  the  simulation  result  (Figure  45a). 
The  frequency  of  pulsing  was  not  changed  and  the  power 
utilized  was  only  negligibly  changed.   The  higher  level 
consumer  (Q5)  was  near  zero  for  the  entire  simulation. 

When  the  model  was  run  with  a  turnover  time  (K17=0.005) 
of  the  top  consumer  (Q5)  longer  than  the  normal  pulsing 
consumer  (Figure  45b),  pulsing  occurred  at  the  normal  fre- 
quency but  the  higher  level  consumer  grew  over  time  until  it 
began  pulsing.   The  period  of  pulsing  became  longer  and  the 
pulse  amplitude  of  the  stored  producer  and  nutrient  storages 
became  greater.   The  normal  pulsing  consumer  (Q3)  remained 
at  a  low  level,  acting  as  a  feeder  to  the  higher  level 
consumer  (Q5).   The  power  utilized  dropped  slightly  to  95.3. 

When  the  turnover  time  of  the  higher  level  consumer 
(Q5)  was  raised  by  another  order  of  magnitude  (K17=  0.0005) 
the  outcome  was  quite  different  (Figure  45c).   The  higher 
level  pulsing  consumer  (Q5)  climbed  toward  an  asymptote 
while  the  stored  production  (Q2)  also  climbed  to  a  steady 
state  value.   the  normal  pulsing  consumer  (Q3)  again  re- 
mained at  a  low  level.   In  this  case  the  nutrients  (Q4) 
became  tied  up  in  the  stored  biomass  (Q2)  the  power  used 
dropped  to  36.1  at  steady  state.   The  percent  power  used  for 
the  entire  simulation  was  48.2. 
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Simulation  of  the  Ring  Model 

The  linear  array  ring  model  was  simulated  with  high 
diffusion  (DK=.l)  between  consumers  in  adjacent  cells 
(Figure  46).   Initially  the  concentration  of  producers  and 
consumers  around  the  ring  was  constant  except  for  a  single 
consumer  at  a  high  level  (Q3  (2,2)  =100.;  lower  left  hand 
corner  of  consumer  matrix).   At  T=50  years  (Figure  46A) ,  the 
consumers  had  pulsed  in  both  directions  around  the  ring  and 
completely  encircled  the  ring  by  T=100  (46B) .   At  T=150 
(46C) ,  the  production  was  beginning  to  spread  around  the 
ring  from  the  lower  left  corner  and  continued  through  T= 
200,  250,  300,  350  (46D-H).   The  consumers  again  began  to 
grow  (T=  350,  46H)  and  spread  around  the  ring  again.   This 
was  followed  by  another  wave  of  production  and  consumption 
(T=  500-750,   46J-0). 

In  runs  with  lower  diffusion  (0.01)  between  consumers, 
the  pulse  wave  traveled  slow  enough  that  the  wave  only  moved 
part  way  around  the  entire  ring  before  the  internal  pulse 
frequency  allowed  the  remainder  of  the  consumers  to  pulse, 
thus  stopping  the  wave.   With  an  even  lower  diffusion  coef- 
ficient of  0.001,  the  wave  moved  3  cells  before  stopping. 
With  a  diffusion  coefficient  of  0.01,  the  wave  moved  10 
cells  before  being  stopped  by  the  natural  internal  pulse 
frequency. 

A  different  pattern  developed  when  the  producers  and 
consumers  in  the  model  were  distributed  in  a  random  pattern 
around  the  ring  (the  individual  cell  concentration  of  pro- 


Figure  46.   Simulation  of  pulse  model  (Figures  14  and  18) 
with  cells  in  a  linear  ring  and  diffusion  between  consumers 
of  each  cell  in  ring  (DK=.l).   For  each  time  unit  (e.g.  A=0) 
density  of  producer  and  consumer  in  the  matrix  is  shown 
along  with  size  class  distribution.   The  time  series  from  A 
to  0  summarizes  the  temporal  pattern  of  totals  in  matrix. 
Initial  conditions  of  consumers  were  set  to  near  zero  except 
for  one  "seed"  consumer  at  lower  left  corner  of  matrix  which 
was  set  to  100. 
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ducers  and  consumers  was  constant  and  the  same  as  the  homo- 
geneous initial  conditions)  and  diffusion  set  to  zero 
(Figure  47).   The  output  was  based  entirely  on  the  random 
field  from  the  initial  conditions.   Each  individual  cell 
model  was  producing  and  consuming  at  the  same  rate  but  there 
was  no  spatial  synchronization  of  the  cells.  The  pattern 
repeated  itself  over  time  (compare  T=50,  47A  with  T=700, 
47N)  . 

When  diffusion  was  set  at  a  high  level  (0.1)  between 
the  consumers,  with  the  same  random  initial  distribution  of 
producers  and  consumers,  the  resulting  pattern  was  quite 
different  (Figure  48).   The  pulsing  consumers  moved  in  a 
wave  around  the  ring  followed  by  a  wave  of  production  (T=50, 
100,  150,  200,  250,  300;  Figure  48A-F)  followed  by  another 
wave  of  consumption  beginning  just  prior  to  T=350  (48G). 
This  was  similar  to  the  simulation  in  Figure  46  that  began 
with  a  homogeneous  initial  distribution  of  producers  and 
consumers  and  had  waves  of  consumption  and  production  around 
the  ring. 

When  the  model  was  run  with  random  distribution  of 
producers  and  consumers  (Figure  49)  and  a  low  value  of 
diffusion  (DK=0.001),  the  spatial  pattern  that  developed  had 
some  properties  of  both  of  the  two  previous  runs.   Because 
speed  of  movement  was  less  with  a  lower  value  of  diffusion, 
a  number  of  focal  points  for  pulse  waves  were  generated 
which  then  run  into  each  other  and  stop.   The  production 
follows  the  pattern  of  consumption  with  multiple  foci. 


Figure  47.   Simulation  of  pulse  model  (Figures  14  and  18) 
with  cells  in  a  linear  ring  but  without  diffusion.   For  each 
time  unit  (e.g.  A=0)  density  of  producer  and  consumer  in  the 
matrix  is  shown  along  with  size  class  distribution.   The 
time  series  from  A  to  0  summarizes  the  temporal  pattern  of 
totals  in  matrix .Initial  conditions  of  producers  and  con- 
sumers were  set  to  random  distribution  around  ring. 
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Figure  43.   Simulation  of  pulse  model  (Figures  14  and  18) 
with  cells  in  a  linear  ring   and  a  high  level  of  diffusion 
between  consumers  of  each  cell  (DK=.l)  and  random  distribu- 
tion of  producers  and  consumers  around  ring.   For  each  time 
unit  (e.g.  A=0)  density  of  producer  and  consumer  in  the 
matrix  is  shown  along  with  size  class  distribution.   The 
time  series  from  A  to  0  summarizes  the  temporal  pattern  of 
totals  in  matrix. 
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Figure  49.   Siniulation  of  pulse  model  (Figures  14  and  18) 
with  cells  in  a  linear  ring  and  a  low  level  ot   diffusion 
between  consumers  of  each  cell  (DK=.0Ol)  and  random  distri- 
bution of  producers  and  consumers  around  ring.   For  each 
time  unit  '(e.g.  A=0)  density  of  producer  and  consumer  in  the 
matrix  is  shown  along  with  size  class  distribution.   The 
time  series  from  A  to  0  summarizes  the  temporal  pattern  of 
totals  in  matrix. 
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Simulation  of  Two  Dimensional  Surface  Models 

The  simplest  simulation  of  the  two  dimensional  pulsing 
del,  with  no  diffusion  and  an  evenly  distributed  energy 
source  (Figure  50),  had  a  time  series  output  identical  to 
the  basic  pulse  model  (Figure  37a).    Even  though  the  model 
was  disaggregated  into  100  cells,  each  of  the  cells  was 
identical.   In  this  run,  each  of  the  cells  was  synchronized 
(by  the  initial  conditions)  and  the  pulsing  was  based  only 
on  the  internal  frequency  of  the  model  (T=250,  50E  and 
T=600,  SOL).   There  was  little  change  in  the  size  distribu- 
tion of  the  producers  and  the  consumers  during  the 
simulation . 

The  influence  of  an  energy  source  that  is  hierarchical- 
ly distributed  from  the  center  of  the  matrix  outward  gen- 
erates a  different  pattern  (Figure  51).   The  production  was 
higher  in  the  center  of  the  matrix  than  at  the  outer  edges. 
In  this  simulation  without  diffusion  there  was  no  edge 
effect.   The  first  pulse  came  at  the  center  of  the  matrix 
(highest  input  energy)  and  then  moved  outward  to  the  edge  in 
a  series  of  pulses.   The  production  and  consumption  then 
continued  to  oscillate.   The  frequency  of  pulsing  in  each 
individual  cell  depended  on  the  intensity  of  the  energy 
input  to  that  cell  (see  also  Figure  40a- f).   The  center 
cells  pulsed  at  a  higher  frequency  than  the  outer  cells  due 
to  differences  in  input  energy.   The  time  series  of  the 


Figure  50.   Simulation  of  pulse  model  (Figure  14)  with 
cells  arranged  in  two  dimensions  (Figure  13)  without 
diffusion  and  with  a  constant  energy  source.   For  each  time 
unit  (e.g.  A=0)  density  of  producer  and  consumer  in  the 
matrix  is  shown  along  with  size  class  distribution.   The 
time  series  from  A  to  0  summarizes  the  temporal  pattern  of 
totals  in  matrix. 
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Figure  51.   Simulation  of  pulse  model  (Figure  14)  with 
cells  arranged  in  two  dimensions  (Figure  18)  .   Energy 
source  hierarchically  is  distributed  from  center  outward  and 
no  diffusion  between  cells.   For  each  time  unit  (e.g.  A=0) 
density  of  producer  and  consumer  in  the  matrix  is  shown 
along  with  size  class  distribution.   The  time  series  from  A 
to  0  summarizes  the  temporal  pattern  of  totals  in  matrix. 
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simulation  had  sharp  peaks  due  to  the  different  frequencies 
of  pulsing  of  the  independent  cells.   The  size  distributions 
of  the  producers  and  consumers  were  based  on  the  input 
energy  and  are  grouped  accordingly.   Without  diffusion,  the 
pattern  formed  was  entirely  dependent  on  the  hierarchical 
pattern  of  the  input  energy. 

The  addition  of  diffusion  between  the  consumers  of  each 
cell  for  the  previous  model  smoothed  out  the  time  series  for 
the  consumers  and  producers  {Figure  52).   A  low  level  of 
diffusion  (DK=0.001)  enabled  the  first  pulsing  cells  (lo- 
cated at  the  center  of  the  matrix)  to  affect  the  neighboring 
cells,  thus  spreading  the  pulse  wave  out  over  the  matrix. 
In  this  simulation  the  size  distribution  of  the  producers 
and  consumers  tended  to  smooth  out  over  time.   The  edge 
effects  were  minimized  in  this  simulation  by  allowing  the 
outer  non-reactive  ring  of  consumer  cells  to  float  at  a 
value  that  was  the  average  of  the  total  consumers  in  the 
matrix . 

Diffusion  between  the  consumers  at  a  low  level  had  a 
much  greater  effect  in  this  two  dimensional  version  of  the 
model  than  in  the  one  dimensional  ring  version  of  the  model. 
When  the  two  dimensional  version  was  run  with  a  random 
energy  source  and  a  low  diffusion  coefficient  (Figure  53, 
DK  =  0.001)  the  effect  was  similar  to  that  seen  in  Figure  52. 
In  this  case,  local  foci  of  high  productivity  (caused  by 
locally  high  values  of  input  energy)  led  to  pulses  that 
spread  over  the  entire  matrix.   This  simulation  was  dif- 


Fiqure  52.   Simulation  of  pulse  model  (Figure  14)  with 
cells  arranged  in  two  dimensions  (Figure  18) .   Energy 
soiree  is  hierarchically  distributed  from  center  outward  and 
diffusion  is  between  consumers  of  each  cell  (DK=.0Ol).   For 
each  time  unit  (e.g.  A=0)  density  of  producer  and  consumer 
tn  thematrix  is  shown  along  with  size  class  -distribution 
The  time  series  from  A  to  0  summarizes  the  temporal  pattern 
of  totals  in  matrix. 
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Figure  53.   Simulation  of  pulse  model  (Uigurs  14)  with 
cells  arranged  in  two  dimensions  (Figure  18) .   Energy 
source  is  randomly  distributed  and  diffusion  is  between 
consumers  of  each  cell  (DK=.0Ol).   For  each  time  unit  (e.g. 
A=0)  density  of  producer  and  consumer  in  the  matrix  is  shown 
along  with  size  class  distribution.   The  time  series  from  A 
to  0  summarizes  the  temporal  pattern  of  totals  m  matrix. 
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ferent  from  the  random  ring  simulations  {Figures  48  and  49) 
in  that  the  input  energy  was  randomly  distributed  while  in 
the  case  of  the  ring  model  the  initial  producer-consumer 
pairs  were  randomly  distributed.   Little  synchronization  of 
the  matrix  occurred  because  the  random  energy  distribution 
caused  locally  high  concentrations  of  producers  every  time 
there  was  a  pulse.   In  the  simulation  of  the  ring  with 
randomly  distributed  producers  and  consumers,  at  a  high 
level  of  diffusion  the  pulse  wave   moved  fast  enough  to 
reset  all  of  the  producers  and  consumers  to  similar  values. 
With  a  low  diffusion  value,  the  wave  traveled  so  slowly  that 
it  did  not  get  around  the  ring,  and  multiple  foci  of  pulsing 
developed . 

Simulation  with  diffusion  between  the  nutrient  com- 
partments of  each  cell  (Q4)  of  the  model  instead  of  to  the 
consumers  (Q3)  can  be  seen  in  Figure  54.   With  a  random 
distribution  of  energy  and  a  high  level  of  diffusion 
{DK=0.1)  the  pulsing  was  almost  totally  uncoupled.   By  the 
end  of  the  run  (T=750,  Figure  540)  there  was  constant  pul- 
sing in  one  cell  or  another,  and  the  overall  level  of  pro- 
ducers as  seen  in  the  time  series  graph  was  fairly  constant. 

The  spatial  configuration  of  the  model  was  also  tested 
with  a  moving  consumer.   This  is  similar  to  the  diffusion 
runs  of  the  model  but  represents  an  active  process  with 
discontinuous  (non-uniform)  movement  of  consumers  from  cell 
to  cell.   The  consumer  was  allowed  to  search  for  the  largest 
producer  to  consume  before  moving.   The  model  was  tested 


Figure  54.   Simulation  of  pulse  model  (Figure  14)  with 
cells  arranged  in  two  dimensions  (Figure  18)  .   Energy 
source  is  randomly  distributed  and  diffusion  between 
nutrient  storages  (Q4)  of  each  cell  is  set  to  high  level 
(DK=.l).   For  each  time  unit  (e.g.  A=0)  density  of  producer 
and  consumer  in  the  matrix  is  shown  along  with  size  class 
distribution.   The  time  series  from  A  to  0  summarizes  the 
temporal  pattern  of  totals  in  matrix. 
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with  a  hierarchical  energy  input  and  a  consumer  search 
length  of  1  cell  (Figure  55)  and  a  search  length  set  to  5 
cells  {Figure  56).   The  simulations  are  quite  different  in 
both  the  spatial  patterns  generated  and  in  the  time  series 
graph  of  the  simulation. 

When  limited  to  a  search  length  of  1  cell,  the  con- 
sumption pattern  moved  like  a  wave  from  left  to  right  across 
the  producers  after  starting  in  the  center.   With  a  longer 
search  length  (Figure  56),  the  consumption  began  in  the 
center  and  spread  out  in  a  circular  pattern  over  the  pro- 
ducers.  There  are  two  of  these  waves  of  consumption  during 
the  time  of  the  simulation  for  the  search  length  of  5.  The 
run  with  a  search  length  of  1  cell  has  slower  consumption 
and  only  moves  across  the  field  once. 

Rain  Forest  Gaps  and  Hierarchies 

Size  Class  Distributions 

Three  different  size  class  distributions  (Figure  57) 
were  generated  from  the  data  set  from  the  radiation  site  at 
El  Verde  to  characterize  the  hierarchical  patterns  in  the 
vegetation.   Figure  57a  represents  the  distribution  of 
plants  by  diameter.  This  can  be  compared  to  the  data  from 
Crow  (1980)  in  Figure  8.   The  distribution  of  plants  by 
crown  diameter  (Figure  57b)  and  by  height  (Figure  57c)  was 
hierarchical.   The  sampling  technique  affected  the  results 
in  the  lowest  size  classes. 


Figure  55.   Simulation  of  the  pulse  model  (Figure  14)  with 
cells  arranged  in  two  dimensions  (Figure  18) .  Moving 
consumer  model  with  search  length  set  to  one  cell,  no 
diffusion  and  hierarchical  energy  distribution.   For  each 
time  unit  (e.g.  A=0)  density  of  producer  and  consumer  in  the 
matrix  is  shown  along  with  size  class  distribution.   The 
time  series  from  A  to  0  summarizes  the  temporal  pattern  of 
totals  in  matrix. 
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Figur-  56.   Simulation  of  the  pulse  model  (Figure  14)  with 
cells  arranged  in  two  dimensions  (Figure  18) .  Moving 
consumer  ^nodel  with  search  length  set  to  five  cells,  no 
diffusion  and  hierarchial  energy  distribution.   For  each 
time  unit  (e.g.  A=0)  density  of  producer  and  consumer  in  the 
matrix  is  shown  along  with  size  class  distribution.   The 
time  series  from  A  to  0  summarizes  the  temporal  pattern  of 
totals  in  matrix. 
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Figure  57.   Size  class  distribution  of  trees  at  El  Verde 
radiation  site  (November  1964) 

(a)  Size  class  distribution  by  diameter 

(b)  Size  class  distribution  by  crown  diameter 

(c)  Size  class  distribution  by  height 
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Gap  Size  Measurements 

The  distribution  of  cecropia  gaps  at  El  Verde  fall  into 
a  hierarchical  distribution  (Figure  58).   Figure  58a  is  the 
size  distribution  of  all  four  of  the  photographic  plots 
combined  and  Figure  58b  shows  the  distributions  of  the 
individual  plots.   The  percentage  of  the  total  area  that  is 
in  the  gap  stage  is  3.79%  (Table  2).   The  values  plotted  in 
Figure  58  are  the  actual  areas  measured  in  square  inches  on 
the  photograph.   The  figure  shows  a  minimum  size  for  the 
gaps  and  a  hierarchical  distribution. 

Comparison  to  Models 

For  each  time  slice  that  the  spatial  simulation  model 
printed  a  spatial  pattern  of  producers  and  consumers,  it 
also  printed  a  graph  of  the  size  distribution  of  the  produ- 
cers and  consumers  (just  to  the  right  of  the  spatial  pat- 
terns).  The  format  of  the  distribution  graph  is  not  the 
same  as  the  size  class  distributions  in  Figure  57  but  the 
size  distributions  do  represent  the  same  class  size  phenom- 
enon.  Depending  on  the  energy  input  conditions  and  dif- 
fusion coefficients  some  of  the  size  distributions  had  simi- 
lar relationships  to  the  natural  distribution  (see  Figures 
51,  53  and  56)  while  others  are  quite  different  (see  Figure 
50).   The  pulsing  in  Figure  50  is  totally  synchronous  while 
the  pulsing  in  Figures  51  and  53  are  more  spatially 
asynchronous. 


ngu.e  58.   Size  distribution  of  Cecropia  gaps  in  tropical 
rainforest  at  El  Verde. 

(a)  Distribution  of  gaps  in  all  five  photographs. 

•I  ,.-^n  of  aaos  in  ^ach  individual  photograph. 

(b)  Distribution  or  gaps  in 


173 


Size    Distribution    of    Cecropia    Gaps 

from  Aorial    Photographs 


SIza   Claaa    (sq.   inchaa) 
1/^X1     Combined    dota 


Size    Distribution    of    Cecropia    Gaps 

from   A«rial    Photographs 


R^    i^^    r.,.^^    PFjri 


0-.02 .02-.04O4-.06O6-.OSOa-.  1  01  0-.  1  21  2-.  1  41  4- .  1  61  6- 

^;.._    <~i *■ 

!7~7\     PICl 


Sizs    Class 

1^^     PIC2  E 


inches) 
PICS 


174 


Table  2.   Area  oE  gaps  digitized  from  photographs  of  Luquillo 
tropical  rain  forest. 


Picture  number 

1 

2 

5 

3 

Total 

Number  of  gaps 

53 

43 

35 

74 

205 

Mean 

0 

.0761 

0 

.0517 

0 

.0455 

0 

.0473 

0.0554 

Std.  Error 

0 

0212 

0 

.0096 

0 

.0043 

0 

.0083 

0.0066 

Minimum 

0 

0093 

0 

.0125 

0 

.0094 

0 

0067 

0.0067 

Ma  X  i  mum 

1 

121 

0 

360 

0 

100 

0 

5525 

1.121 

Area  %  of  total 

5 

38 

2 

.97 

2 

12 

4 

66 

3.735 

*  Means  are  not  significantly  different  p=.005 
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The  gap  size  distribution  was  measured  for  a  set  of 
spatial  simulations  with  input  energy  distributed  hierarchi- 
cally, evenly  and  randomly.   Figure   59a  represents  a  com- 
bined gap  size  distribution  measured  from  a  set  of  three 
different  simulation  runs  using  a  hierarchical  input  energy 
source.   The  gap  size  distribution  is  skewed  to  one  set  of 
large  gaps  and  a  few  smaller  patches.   With  a  random  energy 
source  the  results  (Figure  59b)  resemble  the  size  class 
distribution  of  the  natural  system  (Figure  58)  with  more 
small  patches  and  fewer  large  ones.   With  an  evenly  dis- 
tributed input  energy  source  and  no  diffusion,  the  system 
pulses  in  a  synchronous  manner  that  generates  a  gap  the  size 
of  the  simulation  (100%)  with  each  pulse.   With  diffusion 
present,  the  patch  size  is  dependent  on  the  edge  effect.   If 
the  edge  effect  is  canceled  the  result  is  the  same;  however 
with  a  diffusive  loss  or  gain  along  the  edge,  the  patch  size 
is  reduced  from  100%  due  to  the  uncoupling  of  the  synch- 
ronous pulsing  at  the  edges. 


i^iqure  59.  Size  distribution  of  gaps  in  tropical  rainforest 
pulsing  model  simulation  (Figures  14  and  18)  at  time  =760. 

(a)  Size  class  distribution  from  three  separate  model 
runs  with  hierarchical  energy  distribution. 

(b)  Size  class  distribution  from  three  separate  model 
run  with  random  energy  distribution. 
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CHAPTER  4 
DISCUSSION 

Many  of  the  characteristics  of  ecosystem  function  were 
generated  by  the  simulations  in  this  dissertation.   Energy 
increased  with  growth.  Net  production  alternated  with  pul- 
sing net  consumption.   Hierarchical  patterns  in  space  re- 
sulted  from  oscillations  in  time.   Edge  effects  developed. 
There  were  similarities  with  succession  observed  in  nature. 
Many  characteristics  of  ecosystems  were  generated  by  mini- 
models  that  had  autocatalysis,  recycling,  parallel  pathways 
of  different  order,  spatial  intercell  exchanges  and  hier- 
archical distribution  of  time  constants.   In  other  words, 
simple  models  emulated  many  features  of  more  complex 
ecosystems. 

The  spatial  model  in  this  dissertation  differed  from 
many  previous  spatial  ecosystem  models  that  used  individual 
species  growing  and  interacting  together  (Botkin,  Janak  and 
Wallis  1972,  Phipps  1979,  Doyle  1982).   This  model  was  a 
unit  ecosystem  model  that  combined  all  of  the  species  into 
compartmentalized  production,  consumption  and  nutrient 
storages.   This  simplified  the  model  but  kept  many  of  the 
ecosystem  characteristics. 
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Maximum  Power  Considerations 

The  class  of  models  studied  here  duplicate  real  systems 
by  reinforcing  pathways  that  process  more  power.   The  feed- 
backs simulate  useful  power  processing.   These  models  link 
kinetics  and  energetics  in  ways  observed  in  nature. 

Power  and  Feedback  With  Paths  of  Higher  Order 

Systems  that  generate  higher  order  pathways  to  cap- 
ture varying  energy  flows  may  offer  a  competitive  advantage. 
The  maximum  power  implication  is  that  as  systems  develop 
feedbacks  (higher  order  pathways)  they  can  extract  more 
energy  from  the  source.   Lotka  (1922)  stated  that  as  long  as 
there  was  untapped  available  energy,  systems  were  capable  of 
growth  when  rates  of  flow  increased  through  the  system. 
Odum  (1982  and  1983)  added  that  as  systems  mature  they  feed 
back  energy  which  amplifies  other  pathways  and  maximizes 
power.   The  multiple  pathway  configuration  shown  in  the 
three  path  model  provides  a  possible  mechanism  for  this  to 
occur.   In  the  three  path  model  simulations  (Figures  19  - 
27)  the  linear  pathway  had  a  fixed  efficiency  while  the 
autocatalytic  and  quadratic  pathways  had  variable  effi- 
ciencies (see  Figure  20  and  21)  depending  on  the  input 
power  . 

The  development  of  multiple  pathways  in  a  system  is 
incurred  at  some  energy  cost  to  the  system.   The  energy 
costs  associated  with  developing  and  maintaining  the  non- 
linear pathways  must  be  competitive  to  survive.   For  systems 
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with  small  storages  (i.e.  fast  turnover  times),  the  quad- 
ratic pathway  can  be  non-functional  (Figure  21,  pathway  J3). 
Because  non-linear  pathway  flows  are  a  function  of  both  the 
energy  source  and  the  storage,  there  are  conditions  when  the 
pathway  has  a  threshold  for  operation  (Figure  21,  pathway  J2 
and  J2+J3  and  Figure  27  pathway  J3).   Low  energy  systems  may 
not  have  enough  energy  available  to  allow  development  of 
these  higher  order  pathways. 

Human  systems  may  be  a  good  example  of  how  these  path- 
ways may  operate.   Nomadic,  subsistence  societies  can  be 
considered  as  basically  linear  systems  that  utilize  avail- 
able resources  with  few  or  no  feedbacks.   By  developing 
autocatalytic  feedbacks,  primitive  societies  move  up  to 
developing  societies  building  structures  to  process  more 
energy  (farming,  mining,  transportation  and  manufacturing). 
As  growth  continues,  systems  develop  within  society  that 
have  higher  order  quadratic  feedbacks  to  facilitate  proces- 
sing energy  (communications,  banking  and  finance,  and  infor- 
mation systems).   Because  the  higher  order  pathways  are 
dependent  on  storages  and  energy  flows,  the  structures  may 
not  be  stable  with  reduced  energy. 

For  a  system  pathway  to  utilize  fluctuating  energy 
flows,  it  must  have  enough  structure  to  sustain  the  system 
when  the  non-linear  pathways  are  not  functioning  (at  lower 
energy  levels).   While  the  nonlinear  pathways  were  dependent 
on  the  frequency  and  amplitude  of  input  energy  (Figures  25, 
26  and  27)  the  linear  pathway  had  no  frequency  dependency 
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and  thus  provided  energy  to  the  system  under  all  input 
regimes.   A  system  with  a  combination  of  pathways  then  shows 
greater  stability  under  fluctuating  regimes  and  maximizes 
power  with  increasing  energy  inputs. 

Multiple  pathway  models  have  been  used  to  describe  a 
variety  of  systems.   A  disaster  model  using  multiple  path- 
ways (linear  and  autocatalytic)  has  been  used  to  describe 
earthquakes  and  floods  (Alexander  1978).   Models  of  chemical 
reacting  systems  have  often  used  multiple  pathway  models  to 
describe  the  kinetics  of  the  reactions  ("Brusselator", 
Nicolis  and  Prigogine  1977  and  "Oregonator" ,  Field  and  Noyes 
1974).   "Chaotic  systems"  are  often  modeled  with  multiple 
non-linear  pathways  (Abraham  and  Shaw  1984b). 

Effect  of  Hierarchies  on  Performance 

Hierarchical  subunits  of  a  system  generally  have  in- 
creasing turn-over  times  with  increasing  trophic  levels 
(Allen  and  Starr  1982,  Urban,  O'Neill  and  Shugart  1987). 
The  addition  of  an  extra  consumer  (adding  a  level  to  the 
hierarchy)  of  the  pulse  model  (Figure  15)  must  have  the 
appropriate  turnover  time  to  survive.   If  the  turn-over  time 
was  too  short,  not  enough  energy  was  available  to  that  level 
of  the  hierarchy  to  sustain  it  and  the  added  level  did  not 
survive  (Figure  45a).   If  the  turnover  time  was  too  long, 
the  rate  of  power  use  dropped  and  the  whole  system  collapsed 
(Figure  45c).   The  appropriate  size  consumer  modified  the 
output  behavior  of  the  model  (pulsing  with  a  longer  period), 
but  the  system  was  stable  and  utilized  slightly  more  power. 
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The  highest  level  of  the  hierarchy  in  this  model  determined 
the  frequency  and  scale  of  pulsing.   Therefore  there  are 
optimum  turnover  times  for  maximum  performance. 

Conversely,  as  input  power  increases,  a  higher  level  of 
consumers  may  be  supported.   This  was  seen  in  the  parallel 
production-consumption  model  (Figure  32d)  and  the  pulse 
model  (Figure  40c). 

Power  Used  as  a  Function  of  Input  Power 

The  general  trend  for  all  of  the  models  tested  here  was 
that  as  the  input  power  increased,  the  percent  of  input 
power  that  is  utilized  increased.   This  occurred  in  the 
three  path  model  (Figures  21,  22,  23),  the  parallel 
production-consumption  model  (Figures  30e,  32e,  and  36),  the 
pulse  model  (Figure  40f)  and  the  spatial  models.   This 
appears  to  be  a  function  of  the  non-linear  pathways  that 
feed  energy  back  to  increase  the  efficiency  with  increasing 
available  energy.   Individual  simulations  of  these  models 
with  only  linear  pathways  did  not  show  this  behavior. 

Threshold  for  Stable  Feedbacks  and  Pulsing 

The  pulse  model  exhibited  a  double  threshold  phenom- 
enon. At  low  power  inputs  the  model  did  not  pulse  and  at 
high  power  inputs  the  model  did  not  pulse  (Figures  37  and 
40).  In  the  middle  power  range,  the  model  pulsed  and  the 
pulse  frequency  was  a  function  of  the  input  power.   Local- 
ized maxima  of  power  utilization  may  occur  in  the  pulsing 
range  due  to  synchronization  of  inputs  with  natural  internal 
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frequencies  (Richardson  and  Odura,  1981).   This  double 
threshold  behavior  has  also  been  shown  in  a  wide  variety  of 
prey-predator  model  configurations  (Kuno  1987).   Oscillating 
chemical  reactions  exhibit  this  multiple  output  state  be- 
havior (Field  1985). 

At  low  power  levels,  the  pulsing  model  supported  a 
constant  low  amount  of  consumers  (dependent  on  the  linear 
pathway)  while  at  high  power  levels  the  consumer  was  at  a 
constant  higher  level  (sustained  by  both  the  linear  and 
quadratic  pathways)  with  the  producer  at  a  low  level.   This 
was  also  the  case  in  the  chemical  reactions  and  prey- 
predator  models  described  above.   Models  with  this  behavior 
may  describe  a  variety  of  ecosystems  that  show  various 
levels  of  producers  and  consumers.   A  grassland  ecosystem 
such  as  the  Serengeti  (McNaughton  1985)  may  be  an  example  of 
low  levels  of  producers  supporting  high  levels  of  consumers. 

A  similar  dependence  of  the  highest  trophic  level  on 
the  input  energy  was  also  exhibited  with  the  parallel 
production-consumption  model  (Figure  32)  although  this  model 
did  not  pulse.   It  should  be  noted  with  this  model  that  the 
consumer  level  increased  and  the  'climax'  producer  did  not. 

The  pulsing  model  did  not  pulse  when  the  consumer  quad- 
ratic pathway  (Figure  38)  was  removed,  the  consumer  built  up 
to  a  steady  state,  and  the  percent  power  used  declined. 
There  was  a  lot  of  structure  in. the  higher  level  of  the 
hierarchy  but  the  system  was  not  effective  at  using  the 
extra  power  that  was  available.   Competitively,  a  system 


184 


with  this  structure  may  be  at  a  disadvantage  and  could  be 
eliminated  through  consumption  by  a  higher  level  of  the 
hierarchy  or  competition  by  other  systems  at  the  same  level 
of  the  hierarchy. 

If  a  system  was  not  materials  conservative  (feedbacks 
from  the  consumer  to  the  nutrient  storage  cut  off  or  di- 
verted, Figure  39)  then  the  system  ceased  pulsing  and  ran 
down.  The  system  had  no  feedback  pathways  and  so  did  not 
capture  all  of  the  available  energy. 

Implications  for  Succession 

Role  of  Individual  Units 

Early  successional  producers  can  be  thought  of  as  pre- 
paring the  way  for  succession  to  occur.   Although  early 
successional  species  may  have  other  roles,  in  the  parallel 
production-consumption  model  (Figure  13)  they  can  be  seen  as 
providing  an  energy  source  to  the  consumer  level  of  the 
model  as  the  rest  of  the  system  builds  up.   VJhen  the  early 
successional  species  was  at  a  low  level,  the  consumer  level 
(Q4)  remained  low  (Figure  35).   This  low  consumer  level  did 
not  feed  back  enough  to  the  producers  to  stabilize  the 
system  and  the  system  crashed.   As  the  early  successional 
species  (Q3)  reached  a  threshold  initial  condition,  suf- 
ficient structure  was  built  and  the  system  progressed  to  a 
steady  state. 
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If    the   consumer    level   in   a   successional   system    is   too 
low  then  the  system  may  not  be  stable.      In  the  parallel 
production-consumption   model,    the  consumer   provided  a   feed- 
back on   the   producers    through    the    input   production   multi- 
plier  and    through   consumption   on    the   producers.      When    the 
consumer  was  at  a  level  that  was  too  low,  succession  as 
depicted    by    the    model    (Figure    33   and    36)    did   not   begin.      At 
some   initial   threshold   level  of  consumers,    the   model  pro- 
ceeded   through    a    successional    sequence. 

In  developing   management  plans   for    revegetating    sites 
disturbed   by   mining,    intensive   agriculture   or    natural   dis- 
turbances,   it    is    imperative    that   careful   attention    be   paid 
to  the  whole  structure  of  the  ecosystem  that  is  being  re- 
built.     Without   the  proper   mix  of   early,    middle,    and    late 
successional  producers   along    with   a   set  of  consumers   that 
match   the  producers,    the   restablishraent  of   a   natural   suc- 
cessional   sequence    may   be    retarded   or    destroyed. 

Succession   and   Pulsing 

The  role  of  pulsing   in  succession  may  be  that  in  some 
systems    it    is   necessary   to   have    the   pulsed    recycle    to   main- 
tain  energy   flows  near    maximum    levels.      Several  cases   of   the 
pulsing    model    (Figures    38    and    39)    showed    that    when    recycling 
was   disturbed   power    use   dropped.      Certain    types   of   sucession 
may  need  an  alternation  of   production   and   consumption   at   a 
frequency   that   allows   the   maximum   use  of   available  energy. 
Systems   in   which  available   nutrients   become   bound    in    the 
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bioraass  may  benefit  by  the  fast  release  from  a  pulse  of 
consumption  and  recycle. 

Spatial  Pattern  Formation 

Synchronous  vs.  Asynchronous  Systems 

When  a  spatially  organized  system  is  totally  synchron- 
ized (all  subunits  behaving  as  one),  the  system  may  be  like 
a  monoculture  with  little  pattern  formation  other  than  that 
of  the  local  source  inputs.   In  this  state,  pattern  diver- 
sity is  low.   Where  cells  are  not  all  synchronized  with  each 
other,  patterns  can  develop  that  are  dependent  on  the 
asynchronous  nature  of  the  individual  subunits  as  well  as 
the  local  energy  sources. 

When  the  spatial  model  was  simulated  with  all  of  the 
individual  cells  uncoupled  (not  linked  through  intercell 
diffusion  processes)  and  totally  synchronized  (all  cells 
begun  with  the  same  initial  conditions  and  an  even  energy 
distribution),  no  pattern  was  generated  (Figure  50).   The 
level  of  producers  and  consumers  was  the  same  in  each  cell 
at  every  point  in  time. 

Any  variation  in  the  energy  input  over  the  matrix  area 
lead  to  individual  cells  pulsing  at  frequencies  depending  on 
the  energy  level  local  to  that  area  (Figure  51).   Although 
the  pattern  was  quite  different  from  the  synchronized  one, 
the  energy  use  is  the  same  (Table  5  in  Appendix). 
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Coupling  of  Spatial  Units  by  Diffusion  Processes 

In  any  ecosystem,  spatially  distributed  subunits  are 
connected  to  each  other  through  a  variety  of  processes. 
Nutrients  and  seeds  can  be  carried  spatially  by  transport 
from  wind,  water  and  animal  activity.   Predation  by  con- 
sumers tends  to  reorganize  the  vegetation  community  struc- 
ture.  The  degree  to  which  subunits  are  connected  to  one 
another  is  strongly  reflected  in  the  patterns  that  may 
develop. 

Connectivity  between  subunits  tends  to  decrease  the 
asynchronous  behavior  caused  by  local  energy  differences. 
With  a  low  level  of  diffusion  (Dk=.O01)  the  pulsing  behavior 
was  propagated  across  cell  boundaries  (compare  Figure  52 
with  Figure  51).   At  higher  levels  of  diffusion  (not  shown), 
the  effect  was  to  increase  the  synchronous  nature  of  the 
pulsing  across  the  matrix.   Energy  use  with  various  levels 
of  diffusion  did  not  change  appreciably  (Table  5  in 
Appendix) . 

In  a  single  dimension  system  (ring  model)  the  effect  of 
diffusion  was  similar.   At  high  levels  of  diffusion  (Figure 
48)  pulses  were  propagated  around  the  entire  ring,  while  at 
a  lower  level  of  diffusion  the  propagation  was  confined  to 
local  areas  (Figure  49).   The  asynchronous  pulsing  (Figure 
47)  was  thus  organized  into  a  more  synchronized  spatial 
pattern  depending  on  the  degree  of  connection  between  the 
individual  cells. 
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The  level  of  the  hierarchy  in  which  inter-cell  coupling 
takes  place  plays  an  important  part  in  the  development  of 
spatial  patterns.   When  this  coupling  took  place  at  the 
level  of  nutrient  exchange  from  cell  to  cell  (Figure  54), 
the  effect  was  hardly  detectable,  even  at  high  diffusion 
levels. 

Spatial  patterns  generated  are  not  totally  dependent  on 
the  natural  energy  inputs  but  organize  using  those  natural 
energy  regimes.   Spatial  diversity  thus  depends  on  the  the 
landscape  energy  pattern,  the  interactions  between  the  sub- 
units,  the  hierarchy  level  of  the  interaction  and  the 
existing  pattern  of  vegetation. 

Most  of  the  models  in  this  dissertation  used  only 
diffusive  coupling  between  spatial  subunits  of  the  model. 
Many  systems  have  more  complex  interactions  between  subunits 
than  this  simple  linear  coupling.   The  active  transport 
systems  of  biological  systems  are  good  examples  of  the  more 
complex  coupling  that  can  occur  in  living  systems.   The 
moving  consumer  model  represents  a  more  complex  coupling 
between  individual  cell  units. 

Organization  by  Higher  Level  Consumers 

The  role  of  the  consumer  in  these  models  was  very  im- 
portant in  organizing  pattern  formation.   When  the  spatial 
model  was  simulated  with  one  consumer  spread  evenly  over  the 
matrix,  the  result  was  exactly  the  same  as  when  the  model 
was  simulated  with  all  cells  uncoupled  and  a  single  consumer 
in  each  cell  of  the  matrix  (Figure  50).   In  this  case,  the 
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synchronous   organization   of    individual   consumers    (100    total) 
over   the   entire   matrix  of   cells   mimiced    the   effects   of   a 
large   consumer   with   the   same   territory.      The  percent  power 
used  for  each  of  these  simulations  was  the  same   {Table  5  and 
6  in  Appendix) . 

Coupling  of  the  consumers  from  cell  to  cell  by  dif- 
fusion organized   the  consumer   action  over    the    whole    matrix 
depending   on    the    strength   of    that   coupling.      Low    levels   of 
diffusive   coupling   generated    local   areas   of   organization   by 
the   consumers    (Figure    53)    while    strong   coupling   organized 
the   disturbance   over    the   entire    matrix,    (not   shown   but 
similar    to    Figure    50). 

Active  coupling   between   subunits   by  consumers   was   simu- 
lated  using   a   moving   consumer    model    (Figures    55    and    56).      In 
this  case,    very  different  patterns   were    formed   with   a   smal- 
ler  number   of   consumers.      The   action   of   organizing    the 
entire    landscape    (10x10    matrix)    was   achieved   with   fewer 
consumers.      The   energy   use    was   not   significantly   different 
from    the    other    spatial    simulations    (Table    7    in   Appendix). 
The   efficiency  of   active   coupling    may   be    higher    than   passive 
(diffusion)    coupling. 

Organization   at   a   higher    level   tends    to   have   a    larger 
effect    in   generating   patterns.      Some   of    this    may   due    to   a 
type    of    'memory'    generated    in    the    landscape    by   the 
disturbance-succession    sequence   generated   by    these   pulsing 
production   consumption   models.      As   the   system  pulses,    small 
differences   between   individual   cells   generate    further    dis- 
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continuities.   These  small  differences  act  as  information 
storage  for  future  pattern  development. 

Power  Use  and  Edge  Effects 

No  system  exists  in  an  infinite  plane  without  edges. 
Edges  were  manipulated  in  the  spatial  models  to  understand 
their  role  in  pattern  formation.   Some  of  the  simulations 
allowed  consumers  to  diffuse  into  or  out  of  the  spatial 
matrix  at  high  and  low  levels  of  diffusion. 

When  the  consumer  level  on  the  outside  ring  was  kept  at 
a  low  value  (0.0),  the  percent  power  used  decreased  (Table  8 
in  Appendix)  with  increasing  rates  of  diffusion.   If  the 
outside  buffer  had  a  high  value  for  the  consumer  (Q3  equal 
100)  then  just  the  reverse  was  seen.  With  increasing  rates 
of  diffusion  there  was  an  increase  in  the  percent  power 
used.   This  implies  that  consumer  exchange  can  act  as  an 
energy  source  or  a  drain  in  a  system  depending  on  the 
relationship  of  the  system  to  its  surrounding  area  through 
its  edges. 

General  Principles 

The  following  are  some  general  principles  suggested  by 
model  studies,  which  may  be  useful  hypotheses  in  future 
experimental  studies. 

1)   Multiple  pathways  increase  efficiencies  and  enable 
better  use  of  fluctuating  energy  sources.   Multiple  steady 
states  can  result  from  one  basic  configuration.   The 
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kinetics  of  these  pathway  configurations  are  similar  to 
others  studied  by  chaos  theory,  bifurcation  theory  and 
catatastrophe  theory. 

2)  Hierarchical  structure  is  expressed  in  kinetics  as 
increasing  turnover  times  with  increasing  territory.   Path- 
ways of  control  of  production-consumption  systems  must  match 
the  turnover  time  of  the  appropriate  hierarchical  level  in 
order  to  cause  reinforcement. 

3)  In  early  successional  systems  there  may  be  critical 
minimum  stocks  of  producers  and  consumers  for  a  system  to 
grow. 

4)  Similar  maximum  power  processing  may  be  achieved  by 
a  wide  variety  of  spatial  patterns. 

5)  Connectivity  in  systems  has  a  greater  role  in 
pattern  formation  at  higher  levels  of  the  hierarchy. 
Control  of  patterns  and  patchiness  through  consumer  control 
is  highly  dependent  on  the  spatial  connectivity  of  the 
consumers . 

6)  Patch  size  may  be  related  to  the  turnover  time  of 
the  consumer  and  the  spatial  connectivity  of  the  consumers. 

7)  Some  of  the  great  complexity  of  ecosystems  may  be 
simplified  for  human  comprehension  if  varied  mechanisms  can 
be  grouped  according  to  the  basic  kinetics,  energetics  and 
hierarchical  roles  they  perform. 
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Table  3. 

Coefficient  values  for  parallel  production-consumption 
model  in  Figure  13. 


Production  coefficient  for  Ql 
Production  coefficient  for  Q2 
Production  coefficient  for  Q3 
Drain  coefficient  for  Ql 
Drain  coefficient  for  Q2 
Drain  coefficient  for  Q3 
Consumption  coefficient  for  Ql 
Consumption  coefficient  for  Q2 
Consumption  coefficient  for  Q3 
Drain  coefficient  for  Q4 
Intake  coefficient  for  Q4 
Feedback  loss  coefficient  for  Q4 


Kl 

.003 

K2 

.005 

K3 

.007 

Dl 

.1 

D2 

.2 

D3 

.3 

K7 

.006 

K8 

.015 

K9 

.040 

D4 

.08 

K0 

.1 

Fl 

.01 
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Table  4 

Steady  state  values,  coefficients  and  flows  for  pulse  model. 

Sunlight  normalized  to  100 

Available  sunlight  at  ground  level 

Labile  storage  (Primary  producer) 

Stored  Bioraass 

Pulse  consumer 

Nutrients  (Available  carbon) 


J0     100. 

j;       4.0817993 

Ql    1000. 

Q2   10000. 

Q3      50. 

Q4   30000. 

Kl 

.00000417 

K2 

.5 

K3 

.05 

K4 

.45 

K5 

.00005 

K6 

.00045 

K7 

.0000002 

K8 

.0000018 

K9 

.000002 

K10 

.00000417 

Kll 

.0005 

K12 

.05 

K13 

7.833E-7 

Rl 

510.63309 

R2 

500. 

R3 

50. 

R4 

450. 

R5 

.5 

R6 

4.5 

R7 

5. 

R8 

45. 

R9 

50. 

R10 

510.63309 

Rll 

5. 

R12 

2.5 

(Jordan  and  Drewry  1969,  Odum  and  Pigeon  1970,  and  Brown, 
Lugo,  Silander  and  Liegel  1983) 


Figure  60.   Character  set  for  displaying  spatial  graphs  on 
GIGI  computer  terminal  suitable  for  use  with  screen  copy 
onto  printer.   Each  dot  pattern  is  represented  by  the  hexa- 
decimal code  on  the  left  edge  of  each  plot. 

(a)  80  dots 

(b)  40  dots 

(c)  27  dots 

(d)  20  dots 

(e)  16  dots 

(f)  12  dots 

(g)  7  dots 
(h)   3  dots 
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Table  5.   Percent  power  used  as  a  fjnction  o£  input  energy 
sources  and  diffusion  in  different  ecosystem  levels. 

Successional       Steady  state 

(T)  No  Diffusion  ~~~ 

Energy  Distribution  Percent  Power  Used 

Hierarchical  96.5  96.5 

Even  96.5  96.6 

Random  96.6  96.5 


(b)  Diffusion  between  nutrient  (Q4)  tanks 
(Successional  initial  conditions) 

Diffusion  rate  .001       .01       .1 

Energy  Distribution 

Percent  Power  Used 
Hierarchical  95.5      95.5      95.5 

Even  96.5      "96.5       96.5 

Random  95.6       96.6       95.5 

(c)  Diffusion  between  consumer  (Q3)  tanks 
(Successional  initial  conditions) 

Diffusion  rate  .001      .01       .1 

Energy  Distribution 

Percent  Power  Used 

Hierarchical             95.6       96.6  95.6 

Even                    95.5       96.5  96.5 

Random                  96.6       96.5  95.5 


(a)  Model  DSPl .  See  Figure  50  for  example  run. 

(b)  Model  OSPl.  See  Figure  54  for  example  run. 

(c)  Model  DSPIQZ .   See  Figure  53  for  example  run. 
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Table  6.   Percent  power  used  for  various  runs  of  DSPIC 
spatial  model  having  only  one  consumer  equally  distributed 
across  the  entire  production  matrix. 


Successional  state     Steady  state 


High  initial  condition  for  consumer  Q3  (5000). 

96.5 
96.6 
96.4 

Low  initial  condition  for  consumer  Q3  (50). 

Hierarchical     96.5  96.4 

Even  96.5  96.4 

Random  96.4  96.3 


Hierarchical 

96.5 

Even 

96.5 

Random 

96.4 
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Table  7.   Percent  power  used  for  OSPIOO  model  as  a  function 
of  search  length  and  input  energy  type. 


Search 

length 

Successiona 

1 

Steady 

(cells 

tate 

state 

Hierarchical 

distr 

Lbution 

1 

96.4 

(17) 

(a) 

96 

.5 

(31) 

2 

96.6 

(28) 

95 

.6 

(39) 

3 

96.6 

(36) 

96 

6 

(39) 

4 

96.6 

(36) 

96 

6 

(43) 

5 

96.5 

(33) 

(b) 

96 

5 

(45) 

Even  distribution 

1 
2 
3 

4 
5 

Random  distribution 

1 
2 
3 
4 
5 


95.7 

(15) 

crash 

(13) 

96.4 

(30) 

96.5 

(30) 

96.4 

(38) 

96.5 

(34) 

95.4 

(38) 

96.5 

(42) 

96.5 

(32) 

96.5 

(40) 

96.4 

(20) 

95.7 

(32) 

96.6 

(34) 

96.6 

(3  8) 

96.5 

(40) 

96.5 

(44) 

96.5 

(42) 

95.5 

(45) 

96.5 

(41) 

96.5 

(45) 

(n)  indicates  number  of  consumers  at  end  of  simulation! 

(a)  see  Figure  55 

(b)  see  Figure  56 
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Table  8.   Percent  power  used  as  a  function  of  different 
energy  input  sources  and  diffusion  rates.  Edge  effect  model 
with  different  levels  of  consumers  (Q3)  on  outside  (buffer) 
edge  of  spatial  matrix. 


Diffusion  rate 

.001 

.01 

.1 

Value 

of  Q3 

Energy 

Percent  Power 

Used 

on  outside 

Distribution 

edge 

Hierarchical 

96.6 

96.5 

96.1 

0.0 

Even 

96.5 

95.5 

95.9 

Random 

96.6 

96.5 

95.9 

Hierarchical 

96.6 

96.5 

96.7 

50. 

Even 

96.6 

96.7 

96.5 

Random 

96.6 

96.7 

96.6 

Hierarchical 

96.6 

96.7 

96.3 

100. 

Even 

96.6 

96.8 

96.8 

Random 

96.6 

96.3 

96.3 
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PROGRAM  sue 10 
C        SUCGGX 
C        VERS  1  . 1 
C        FEBRUARY  5,  1984 


C 


BYTE  F1LE(16),ESC,DES(40) 
REAL  Ml ,M2,M3,M4,M9 

REAL  K1,K2,K3,K4,K5,K6,K7,K8,K9,K0,L1,L2,L3,J,J0 

DIMENSION  FILNAM(6) ,ly(50,200) 

DATA  PILE/16*0/ 

DATA  DES/40*0/ 

FT1(A,B)=ABS(AINT(A/B)-A/B) 

D(X,Y)={X/y)*ALOG(X/Y) 


WRITE(5,100) 

100  FORMATdX,'       SUCGGM  GENERATES    6    DATAFILES' , 
S/'       BE   SURE   THAT   THEY   DONT   ALREADY    EXIST', 

S/'       WHAT    IS   THE   DATA   FILE   FOR  THIS    MODEL    RUN    7 ' ) 
READ(5,101)(FILE(I),I=1,16) 

101  F0RMAT(16A1) 
C 

C 

C  WRITE(5,1011) 

C1011         FORMATC    WHICH   Q  TO   SAVE    (1,2,3,4,5=%    POW   USED,6=BI0MSS) ' S) 

C  READ(5,1012)IQSAV 

C1012         F0RMAT(I3) 

WRITE(S,1013) 

1013  FORMATC       WHAT    IS   THE   INCREMENT    IN   J07     [  R]     •$) 
READ(5,1014)XINC 

1014  FORMAT{G15.5) 
C 

C 

c        WRITE(5,99) 

c  99     FORMAT( '    HOW  LONG  TO  RUN?  •) 

c        READ(5,98)TIME 

o  98     FORMAT(P6.0) 

TIME=100.  IX. 1 

C        WRITE(5,981) 

C981     FORMATC   DO  YOU  WANT  A  HARDCOPY?  (  1-YES,0-NO)  '  S ) 
C        READ  ( 5,982) ICOPY 
C982     FORMAT (12) 

0PEN(UNIT-1, NAME-FILE, TYPE='OLD' ,FORM='UNFORMArTED' ) 

READ(1)E1,E2,E3,E4,ES,E6,E7,E8,E9,E10,E11,E12,E13,E14,E15 
C        COEFFICIENTS  *****♦*♦♦♦**♦*«****♦*«*♦*»•*♦***»♦«, «*«,«*,»^, 

+NUM,K0,K1,K2,K3,K4,K5,K6,K7,K8,K9,D1,D2,D3,D4,L1,L2,L3,F1 
+,J0,Q1INIT,Q2INIT,Q3INIT,Q4INIT 
CL0SE(UNIT=1) 
C        INITIAL  C0NDITI0NS**»********»***»««**««»****«*«*»»*,,»»,4,»«, 
XJ0INI=J0 
NSLICE=25 
NCNTS=150 
1        CONTINUE 
C 
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C        BIG  OUTER  LOOP 
C 

DO  1062  IQSAV=1,6 
C 
C 

NRUN=0 
2        CONTINUE 

JO=XJOINI+NRUN*XINC 

NRUN=NRUN+1 

T=0 

PERCNT=0 

PRVAII;=0 

PUSED=0 

DUSED=0 

Q9=0 

BIOMSS=0 

M1=0 

H2=0 

M3=0 

M4=0 

M9=0 

P=0 

BMAX=0 

Q1=Q1INIT 

Q2-Q2INIT 

Q3=Q3INIT 

Q4=Q4INIT 

Q1S1ZE=30. 

Q2SIZE=5. 

Q3SIZB=1. 

Q4S1ZE=20. 
C        SET  OUTPUT  VECTOR  AND  FLAG  ***♦****»*♦**♦**♦**♦***♦*****♦•«»« 

DT=.  1 
C        WRITE(5,108)DT 

C  108    FORMAT ( '   TIME  INTERVAL  DT=  ',F5.3) 
C        ISTEP=1/DT  IDr'S  PER  T 

C        IPLOT=0  [PLOTTING  INTERVAL  IN  DT' S 

C        ITCNT=0  ! ITERATION  COUNTER 

C        WRITE(5,1081) 

C1081    FORMATC  WHAT  IS  PLOTTING  INTERVAL  PER  TIME  UNIT  [I]  ■$) 
C        READ(5, 1082)IPLOT 
C1082    F0RMAT(I2) 
C        HRITE(5, 1083) 

C1083    P0RMAT(1X,'  INPUT  VALUES  FOR  SIZES  Q1-24  [ R]  'S) 
C        READ  (5,1084)Q1SIZE,Q2SIZE,Q3SIZE,Q4SIZE 
C1084    PORMAT(4F8.3) 

C        WRITE  GIGI  STARTUP  INFORMATION 
C        CALL  GGON 
C        CALL  GGINIT 
C        CALL  GGERA 
C        CALL  GGAXIS(0,0,767,479) 
C        CALL  GGB0X(7, 0,0, 767, 479) 
C        CALL  GGBOX(7,0,0,767,350) 
C 
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*♦**♦*♦**♦***••**»***•••♦«******************* 


C  XDT=DT/10. 

C  START   OF   LOOP   *********************************************** 

5                 T=T+DT 
C                   1TCNT=ITCNT+1 
C  RATE   EQUATIONS*"*****""*"""**"***"***""***** ** 

J=JO/( 1+L1*Q1*Q4+L2*Q2*Q4+L3*Q3*Q4) 

EUSED=JO-J 

R1=DT*K1*Q1*Q4*J 

R2=DT*K2*Q2*Q4*J 

R3=DT*K3*Q3*Q4*J 

R4=DT*Q1*D1 

R5=DT*Q2*D2 

R6=DT*Q3*D3 

R7=DT*K7*Q1*Q4 

Ra=DT*K8*Q2*Q4 

R9=DT*K9*Q3*Q4 

R0=DT*D4*Q4 
C  LEVEL    EQUATIONS 

Q1=Q1+R1-R4-R7 

Q2=Q2+R2-RS-R8 

Q3=Q3+R3-R6-R9 

Q4=Q4+K0*  (  R7+R8+R9 )  -R0-F1*  (  R1  +  R2+R3  ) 

Q9"R1+R2+R3+K0*(R7+R8+R9) 

DRAIN-(R4+R5+R6+R0)+(  1  .-KO)  *(  R7+Ra+R9) 

BIOMSS=Q1+Q2+Q3+Q4 
C  M1=AMAX1(M1 ,Q1) 

C  M2=AMAX1(M2,Q2) 

C  M3=AMAX1(M3,Q3) 

C  M4=AMAX1(M4,Q4) 

C  M9=A;'IAX1(M9,Q9) 

c 
c 

C  P1=Q1/Q1SIZE 

C  P2=Q2/Q2SIZE 

C  P3=Q3/Q3SIZE 

C  P4=Q4/Q4SIZE 

C  P1-AMAX1(P1 ,1E-5) 

C  P2=AMAX1(P2,1E-5) 

C  P3=AMAX1(P3,1E-5) 

C  PMAX=P1+P2+P3+P4 

C  DIVERS=-(D(P1,PMAX)+D(P2,PMAX)+D(P3,PMAX)+D(P4,PMAX)) 

C 

D  WRITE(5,103)P1,P2,P3,P4,PMAX 

D103  FORMATC 1X,5(2X,G12.5) ) 

C 

C 

C  BMAX=AMAX1(BMAX,BI0MSS) 

P=P+Q9 

PAVAIL=PAVAIL+JO*DT 

PUSEI>=PUSEEH-EUSED*DT 

DUSED=DUSED+DRAIN 

PERCNT= BUSED/ JO 
C 
C  FIND    WHICH   Q   TO   SAVE 
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20000   GOTO( 2 1000, 22000, 23000, 24000, 25000, 26000 )IQSAV 

GOT01101 
21000    IY(NRUN,INT(T*1.5+1))=INT(Q1/.5) 

GOT0 1101 
22000    1Y(NRUN,INT(T*1.5+1) )=INT(Q2/.5) 

GOTO1101 
23000    IY(NRUN,INT(T*1.5+1) )=INT{Q3/.5) 

GOTO1101 
24000    IY(NRUN,INT(T*1.5+1) )=INT(Q4/.5) 

GOT01101 
25000    IY(NRUN,INT(T* 1.5+1) )=INT( PERCNT* 1000 . ) 

GOTO1101 
26000    IY(NRUN,INT(T*1.5+1) )=INT(BI0MSS/. 5) 

GOT01101 
C 

C        SKIP  PLOTTING  IN  THIS  VERSION 
C 
C 

C        GOTO  1101 

C        IF(FT1(T,.1).GE.DT)G0T0  1101 
C        IF(ITCNT.LT.IPLOT)GOT01101 
C        ITCNT=0 
C        IX=T*7. 
C        IY=Q1 

C        CALL  GGPLT(6,IX,IY,1)    [YELLOW  FOR  CLIMAX  SPECIES 
C        IY=Q2 

C        CALL  GGPLT( 1,IX,IY,1)    IBLUE  FOR  TRANSITIONAL  SPECIES 
C        IY=Q3 

C        CALL  GGPLT(2,IX,IY,1)    !  RED  FOR  WEEDS 
C        IY=Q4 

C        CALL  GGPLT(3,IX,IY,1)    ! MAGENTA  FOR  CONSUMERS 
C        IY=Q9/DT 

C        CALL  GGPLT(4,IX,IY,1)    1  GREEN  FOR  PRODUCTIVITY 
C        IY=BIOMSS 

C        CALL  GGPLT(5,IX,IY,1)    ICYAN  FOR  BIOMASS 
C 
C 

C        I Y= BUSED* 100/ JO  'SCALE  EUSED  TO  0-100 

C        CALL  GGPLT(7,IX,IY+350,1)  IPLOT  POWER  USED   WHITE 
C        IY=(DRAIN/DT)*100/J0       'SCALE  DRAIN  TO  0-100 
C        CALL  GGPLT(5,IX,IY+350,1)  ICYAN  FOR  DRAINS 
C        IY=DIVERS*50 
C        CALL  GGPLT(2,IX,IY+350,1) 
1101     IF(T.LT.TIME)GOTO  5 

WRITE( 5 , 1 10 1 1 ) JO ,Q1 ,Q2 ,Q3 ,Q4 , PERCNT* 100 . ,BIOMSS 
11011    FORMAT( 1X,7( 1X,F10.4) ) 

IP ( NRUN . LT . NSLICE ) G0T02 

CALL  ASSIGN(2, 'SUCMANY' ,6) 

ENCODE (40, 25001, DBS )IQSAV, FILE 
25001    FORMATdX,'  TANK  Q", II,'  FOR  DATA  PILE  ',16A1) 

WRITE(2) (DES(I) ,1=1,40) , NSLICE, NCNTS, 
+   ((IY(JCNT,KCNT) ,KCNT=1 , NCNTS) , JCNT=1 , NSLICE) 

CL0SE(UNIT=2) 

WRITE(5,1061)IQSAV 
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1061  FORMATC       END   OF   RUN   #    ',14) 

1062  CONTINUE 

C  ESC=27 

C  IF    (lCOPY.EQ.0)GOTO1102 

C  WRITEO, 11021) 

C11021       FORMATC     S(H)") 

C  CALL  GGERA 

C1 102  CALL  GGOFF 

C  WRITE(5,1009)ESC,ESC,ESC 

C1009  FORMAT{ 1X,A1 , 'Prtml' , A1 , '«' , A1 , ' [H'// 

C  &       '       INITIAL  VALUES   OF  VARIABLES') 

C  WRITE(5,114) 

C  WRITE (5,1 12) Q1INIT,Q2INIT,Q3INIT,Q4INIT,0.,0. 

C  WRITE(5,111) 

cm  FORMATC         MAXIMUM  VALUES   OF  VARIABLES       ') 

C  WRITE(5,114) 

C    114  F0RMAT(6X,'Q1' ,7X, 'Q2' ,7X, 'Q3' ,7X,'Q4' ,6X, 'PROD' ,5X,'BI0MASS' ) 

C  HRITE( 5, 1 12)M1 ,M2 ,M3 ,M4 ,M9/DT,BMAX 

C    112  FORMAT( '     ' ,6F9.3) 

C  WRITE(5,113) 

C    113  FORMATC          FINAL   VALUES   OF  VARIABLES') 

C  WRITE(5,114) 

C  WRITE(5,112)Q1,Q2,Q3,Q4,Q9/DT,BIOMSS 

C  WRITE(5,115)PAVAIL,P,PUSED,(PUSED/PAVAIL)*100,DUSED 

C    115  FORMAT( 

C  Six,'    ENERGY   AVAILABLE   =          ',F12.4/ 

C  Six,'    TOTAL   PRODUCTIVITY   =    ',P12.4/ 

C  Six,"    ENERGY   USED   =                       ',P12.4,'    PERCENT   USED   =       ',F12.4/ 

C  SIX,'     ENERGY   DRAINED  =               ',F12.4) 

C  WRITE(5,116) 

C    116  FORMAT(6X,'K0' ,8X, 'K1' ,8X, •K2' ,8X,'K3' ,8X, 'K7' ,8X, 'K8' ,8X, 

C  +'K9' ,8X, 'F1' ) 

C  WRITE(5,117)K0,K1 ,K2 ,K3 ,K7 ,K8,K9,F1 

C    117  FORMAT ( '     ' ,8( 2X ,F3 .4) ) 

C  WRITE(5,118) 

C    118  P0RMAT(6X, 'D1' ,8X,'D2' ,8X, 'D3' ,8X,'D4' ,8X,'L1' ,8X,'L2' ,8X, 

C  +'L3' ,8X,' JO' ) 

C  HKITE(5,117)D1 ,D2 ,D3 ,D4 ,L1 ,L2 ,L3 ,J0 

C  WRITE(5,119)DT,TIME 

C  119  FORMAT ( '   DT  THIS  RUN  =  ',F6.4,'  TOTAL  T=  ',F6.2) 

C  WRITE(5,120) (FILE(I) ,1=1 ,16) 

C  120  FORMATC   DATA  FILE  DESIGNATION  FOR  THIS  RUN  ',16A1) 

C  WRITE(5, 121)Q1SIZE,Q2SIZE,Q3SIZE,Q4S1ZE 

C121  FORMATC 3X,'Q1 SIZE  Q2SIZE  Q3SIZE  Q4SIZE'/ 1X,4F7 . 1) 

C  IF(ICOPY.EQ.0)GOTO1201 

C  CALL  GGON 

C  WRITE(3,11021) 

C  CALL  GGOFF 

1201  END 
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1  REM  THREPftTH    MODEL  VERSION   7/1/34 

2  REM      WITHINPUT   COEFFICIENTS   CHANGED 

3  REM  NAME=>    TPM0D7.BAS 

90  PRINT    "WHAT    IS   VALUE   FOR   ENERGY    INPUT    (J0-J1)" 

91  INPUT   J1 

100    REM  AND   PRINT    PNVS(C,X,Y)    FOR  PLOTTING   VECTORS.    C=C0L0R(1-7) 
110    PRINT    "DO   YOU   WANT   GRAPHICS    ON"    •    INPUT   Q$ 
120    IF   Q$<>"Y"    GO   TO    150 

130    PRINT   CHR5(27)+"PpS(E)H(R,I(G) ,P1 ,NO,AO ,S0 ) S( A[ 0 ,479] [767,0] )" 
140    DEF   FNV$(C,X,Y)="W(I"+STR$(C)+")V("+STRS(X)+","+STR$(Y)+"]" 
150    DEF   PNP5(C,X,Y)="W(I"+STR$(C)+")P["+STRS(X)+","+STR5( Y)+"1V[]" 
160    DEF   FNT$(C,N,A$)="W(I"+STR$(C)+")T(S"+STR$(N)+") '"+A$+"'" 
170    DEF   FNB$(C,X,Y,X1,Y1)=FNPS(C,X,Y)+FNV$(C,X,Y1)+FNV$(C,X1 ,Y1) 

+PNV$(C,X1,Y)+PNV$(C,X,Y) 
180    A$=CLK$ 

190  A9=TTYSET(2S5,132) 
200  N9=-1.8 
210  T9=1 
220  H=0 
230  5=100 
240  K1-.5 

250  K2=1.00000E-03 
260  K3=1 .OOOOOE-06 
270  K4=.2 

275  K6=K1 

276  K7=10*K2 

277  K8=10*K3 
280    T=0 

290   A$="Threepath  Model" 

300    PRINT    FNPS(7,626,475) ;FNT$(7,1 ,A5) 

310    PRINT    FNB$( 7,620,456,767,479) 

330   PRINT   PNB$(7,0,0,767,479) 

340    PRINT    PNPS(7,0,270) ;FNV$(7,767,270) 

350    PRINT    FNP$(7,0,380);FNV$(7,767,380) 

360    PRINT    FNP$(7,0,244) ;FNVS(7,767,244) 

370    J0=J1/2+COS(W*T/57.2958)*J1/2 

380    J9=J0/(  1+K6+K7*!J+K8*Q*Q) 

390    P1=P1+J0-J9 

400    J7=J7+J0 

410    R1=K1*J9 

420    R2=K2*Q*J9 

430    R3=K3*Q*Q*J9 

440    R4=K4*2 

450    Q9=T9*(R1+R2+R3-R4) 

460    Q=Q+Q9 

470    X0=R1+R2-I-R3 

480    X1=100*R1/X0 

490    X2=100*R2/X0 

500    X3=100*R3/X0 

510    IF   Q$="Y"    THEN    580 

520    PRINT    "J0=";J0,"J9=";J9,"T=";T 

530    PRINT    "S=",-Q,"Q9=";Q9 

540   PRINT   "R1=";R1,"R2=";R2,"R3=";R3,"R4=";R4 

550    PRINT    "X1=";X1,"X2=";X2,"X3=";X3 
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560   PRINT 

570    IF   QS<>"y"    THEN   640 

580    PRINT    FNP$(1,T*2,X1+275),FNP$(2,T«2,X2+275) 

,FNP$(3,T*2,X3+275),PNP$(4,T*2,5+J9/1) 
590    PRINT    FtIP$(5,T*2, JO/100+5) 
600    PRINT    FNP$(7,T*2,Q/1000+385) 
610    PRINT    FNP5(4,600,270) !"T(S1)'JR="!J9;"'" 
620    PRINT    FNP$(6,T*2,170+( 100*J9/JO) ) 
630    Q7=27+Q  ®    REM  TOTAL   Q   TO   GET    AVERAGE 
640    T=T+T9  ®   IP   T<360    THEN    370 
650   PRINT   FNP$(4, 600, 270) ;"T(S1)' 
660   J5=P1/J7*100 
670    PRINT    FNP$(7,0,270);"T(S1)'POW  USED=";P1 

,-"POW   AVAIL=";J7;  "PERCENT   OSED="  ;  J5  ;"AVE   2="!Q7/T;" 
680    INPUT   X 
690    PRINT   CHR$(27)+"»" 

700    PRINT    P1/J7; "FRACTION   OF  TOTAL   POWER  USED" 
710    END 


c 


203 


c 

c 

C        GIGI  GRAPHICS  SUBROUTINE  PACKAGE 

C        WRITTEN  BY  JOHN  R.  RICHARDSON 

C 

C        SEPTEMBER  1982 

C 

C  VERS  11:   ALL  UPDATES  AND  CURRENT  TO  SEPTEMBER  1982 

C 

C  VERS  12:  FEBRUARY  28  1984  ADDITIONS 

C  ADDED  GGPLOT  (CALL  TO  GGPLT) 

C  ADDED  GGDMP   ( HARDCOPY  DUMP) 

C  ADDED  GGVERS 

c 

C        ALL  I/O  IS  TO  LOGICAL  UNIT  3 
C 

C        THE  NORMAL  CALLING  SEQUENCE  TO  SET  UP  THE  GIGI  WOULD  BE         C 
C        AS  FOLLOWS:  ^ 

C 

C        CALL  GGON  'TURNS  GRAPHICS  ON  C 

C        CALL  GGINIT  ! SENDS  NORMAL  INITIALIZATION    C 

C        CALL  GGERA  'ERASE  THE  SCREEN  C 

C        CALL  GGAXIS(0, 0,767, 479)         !SETS  NORMAL  AXIS  WITH  ORIGIN   C 
'^  !  AT  THE  LOWER  LEFT  CORNER      C 

C  TASKBUILDING  USING  THE  GIGI  ROUTINES  C 

C  RUN  THE  TASK  BUILDER  (TKB  <CR>)  C 

C  TKB>MYPRDG=MYPROG,LB: [ 1,1]GGLIB/LB  <CR>  C 

C  TKB>/  J, 

C  ENTER  OPTIONS  2 

C  TKB>ASG=TTn:3  IWHERE  n  EQUALS  GIGI  TERMINAL  NUMBER    C 

C  TKB>//                   !   COULD  USE  TI:  INSTEAD  C 

C 
C 
C 

SUBROUTINE  GGVERS  ( IVERS  ) 
C        CALL  TO  THIS  WILL  GIVE  THE  CURRENT  VERSION  OF  THE  GGLIB 

IVERS=12 

RETURN 

END 
C 
C 
C 

SUBROUTINE  GGON 
C        THIS  WILL  SEND  THE  ESC  Pp  SEQUENCE  TO  THE  GIGI  TO  ENABLE  THE 
C        GRAPHICS 

BYTE  ESC 

ESC=27 

WRITE(3,100)ESC 
100      FORMAT( '  +  ' ,1A1, 'Pp' ) 

RETURN 

END 
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C 

c 
c 

SUBROOTINE  GGOFF 
C        THIS  WILL  SEND  THE  ESC  »  NEEDED 
C         TO  TURN  OFF  THE  GIGI  GRAPHICS  MODE 

BYTE  ESC 

ESC=27 

WR1TE(3,100)ESC 
100      FORMATt '+' ,1A1, ■«' ) 

RETURN 

END 
C 
C 
C 

SUBROUTINE  GGERA 
C        ROUTINE  TO  PERFORM  SCREEN  ERASE 

WRITE  (3,100) 
100      PORMAK '  +  ' , 'S(E)  '  ) 

RETURN 

END 
C 
C 
C 

SUBROUTINE  GGDMP 
C  ROUTINE  TO   PERFORM   SCREEN   DUMP   TO   LA34/LA10D    PRINTER 

WRITE(3, 100) 
100  FORMAT('+' , 'S(H) ' ) 

RETURN 

END 
C 
C 
C 

SUBROUTINE   GGINIT 
C  ROUTINE   TO    INITIALIZE   THE  GIGI 

WRITE(3,100) 
100  FORMAT ( '+' , •W(R,I4,P1 ,NO,SO,AO) ' ) 

RETURN 

END 
C 
C 
C 

SUBROUTINE   GGAXIS( IX.IY ,IFX, IFY) 
C 

C  ROUTINE  TO   INITIALIZE   THE   AXIS   OF  THE   GIGI 

C  WHERE   IX  =   LOWER  LEFT   CORNER  X   VALUE 

C  lY  =    LOWER  LEFT   CORNER  Y  VALUE 

C  IPX=    UPPER   RIGHT   CORNER  X   VALUE 

C  IFY=    UPPER   RIGHT   CORNER  Y  VALUE 

WRITE(3,100)IX,IFY,IFX,IY 
100  FORMAT ( '+' , 'S(A[" ,15, ■ , ' ,15, ']  [' ,15,  ■ , ' ,15, •]  ) ' ) 

RETURN 

END 
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C 

c 

SUBROUTINE  GGPLT(COLOR,IX,iy ,IPLAG) 
C   SUBROUTINE  TO  POSITION  GRAPHICS  CURSOR  ON  THE  GIGI 
C     SET  IFLAG  TO  NUMBER  >  0  TO  PLOT  POINT  AND  TO  0 
C     TO  MOVE  CURSOR  TO  POINT  WITHOUT  PLOTTING  POINT 
C  COLOR  =  BYTE  variable  0-7  for  color 

C  IX  =  Integer  value  of  X 

C  lY  =  Integer  value  of  Y 

C  IPT  =  Integer  flag  >1  plot  a  point 

BYTE  COLOR 

IF( IFLAG. GT.O) GOTO 10 

WRITE(3, 100)COLOR,IX,IY 

100  FORMAT('+  H(I' ,I1,')P1' ,14,', ',14,'] ') 
GOTO20 

10       WRITE(3,101)COLOR,IX,IY 

101  FORMAT( •+  W(I' ,11, ' )P[ ' ,14, ', ' ,14, "IVtl ■ ) 
20       CONTINUE 

RETURN 

END 
C 
C 
C 

SUBROUTINE  GGPL0T(C0LOR, IX, lY, IFLAG) 
C  ROUTINE   TO   ALLOW   FOR  VARIATION    IN   SPELLING   OF   GGPLT    ROUTINE 

C  ADDED    IN  VERS    12 

BYTE   COLOR 

CALL  GGPLT (COLOR, IX, I Y, I FLAG) 

RETURN 

END 
C 
C 
C 

SUBROUTINE  GGVEC (COLOR, IX, lY) 
C   SUBROUTINE  TO  DRAW  A  VECTOR  ON  THE  GIGI  FROM  ITS  PRESENT  POSITION 
C     TO  THE  IX, I Y  POSITION  IN  THE  PARAMETER  LIST.   USE  GGPLT  FOR 
C     INITIAL  COORDINATES  IF  NEEDED. 
C  COLOR  =  BYTE  variable  0-7  for  color 

C  IX  =  Integer  value  of  X 

C  lY  =  Integer  value  of  Y 

BYTE  COLOR 

WRITE (3,100) COLOR, IX , lY 
100  FORMAT( '+  W(I' ,11, ')V[ ' ,14, ' , ' ,14, '] ' ) 

RETURN 

END 
C 
C 
C 

SUBROUTINE  GGBOX ( COLOR, IX , I Y , IX 1 , IY1 ) 
C   SUBROUTINE  TO  DRAW  A  BOX  ON  THE  GIGI  GIVEN  THE  OPPOSITE  COORDINATE 
C     PAIRS  FOR  THE  RECTANGLE. 

C  COLOR  =  BYTE  variable  0-7  for  color 

C  IXO  =  Integer  value  of  X 

0  lYO  =  Integer  value  of  Y 
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C  1X1  =  Integer  value  of  X  opposite 

C  lYI  =  Integer  value  of  Y  opposite 

C 

C        IF  THE  FILL  IS  TURNED  ON  THE  BOX  WILL  BE  FILLED  AUTOMATICALLY 

C 

BYTE  COLOR 

CALL  GGPLT(C0L0R,IX,IY,1) 

CALL  GGVEC(C0L0R,IX1,IY) 

CALL  GGVEC(C0L0R,IX1,IY1) 

CALL  GGVEC(C0L0R,IX,IY1) 

CALL  GGVEC(COLOR,IX,IY) 

RETURN 

END 
C 
C 

c 

SUBROUTINE  GGCIRC(COLOR,IX,IY,IRAD) 
C   SUBROUTINE  TO  DRAW  A  CIRCLE  AT  POINT  IX, I Y  WITH  A  RADIUS  OF  IRAD 
C  COLOR  =  BYTE  variable  0-7  for  color 

C  IX  =  Integer  value  of  X 

C  lY  =  Integer  value  of  Y 

C  IRAD  =  Integer  radius  of  circle 

C 

C   IF  THE  FILL  IS  TURNED  ON  THE  CIRCLE  WILL  BE  PILLED  AUTOMATICALLY 
C 

BYTE  COLOR 

CALL  GGPLT(COLOR,IX,IY,0) 

WRITE(3,100)COLOR,IX,IY+IRAD 
100      FORMAT('+  W(I' ,I1,')C[' ,14, ' , • ,14, 'J ' ) 

CALL  GGPLT(C0LOR,IX,IY,1)  [LEAVE  CURSOR  AT  CENTER 

RETURN 

END 
C 
C 
C 

SUBROUTINE  GGTEXT( COLOR, IX , I Y ,TEXT , ISI ZE , ITILT ) 
C  SUBROUTINE  TO  WRITE  TEXT  AT  IX, I Y  ON  SCREEN 
C       Writes  text  at  ix,iy  with  size  and  rotation  of 
C        characters  given 

C  COLOR  =  BYTE  variable  0-7  for  color 

C  IX  =  Integer  value  of  X 

C  lY  =  Integer  value  of  Y 

C  TEXT  =  BYTE  array  containing  80  char  or  less 

C  ISIZE  =  Integer  value  for  text  size  0-8 

C  IROT  =  Integer  value  of  degress  of  rotation  for 

C  line  of  text  (multiple  of  45) 

C 

C        THIS  SUBROUTINE  CALLS  LENGTH  TO  DETERMINE 
C         THE  LENGTH  OP  THE  STRING 
C 

BYTE  COLOR, TEXT( 1) 

N=0 

CALL  GGPLT(COLOR,IX,IY,0) 

CALL  LENGTH(TEXT,N) 
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WRITE    (3,100)ITILT,ISIZE,ITILT,(TEXT(I) ,I=1,N) 
100  FORMAT('+   T(D' ,I4,')(S' ,12, ')(D' ,14,') ' ,1H',<N>A1,1H') 

RETURN 

END 
C 
C 

c 

SUBROUTINE   LENGTH ( TEXT , N ) 

BYTE  TEXT(80) 

IFLAG=0 

DO   20    1=80,1,-1 

IP(TEXT(I).GT.32)IFLAG=1 

N=I 

IF    (IPLAG.EQ.0)GOTO20 

G0TO99 
20  CONTINUE 

99  RETURN 
END 

C 
C 

c 

SUBROUTINE  GGPILL( IFLAG) 
C        SUBROUTINE  TO  TURN  ON/OFF  COLOR  FILL  CHARACTERISTIC 
C  IFLAG=0  NOPILL,  IFLAG=1  FILL 

C 

WRITEO,  100)  IFLAG 

100  FORMAT( '+' , 'W(S' ,11, ' ) ' ) 
RETURN 

END 
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PROGRAM  PLOTZ 

c 
c 

C  **6/27/83 

C        CHANGED  AXIS  ROUTINE  FOR  THREECORNERED  ORIGIN** 

C 

C        VERSION  1.6 

C        WRITTEN  BY  JOHN  RICHARDSON 

C        APRIL  27,  1983 

C 

C        SURFACE  PLOTTING  PROGRAM 

C 

DIMENSION  IY(50,200) ,IOUT(200) ,IOUTY(200) 

DIMENSION  IX(200),MASK(a00) 

BYTE  FNAM(20) ,DES(40) ,G0N(3) ,G0FF(2) , COLOR 

BYTE  BLACK , BLUE , RED , MAGENT , GREEN , CYAN , YELLOW , WHITE , ESC 

INTEGER  DELTAXjDELTAY 

COMMON  /I AREA/MASK 

DATA  FNAM/15*0,' . ' , 'D' , ' A' , 'T' ,0/ 

DATA   MASK/800*0/,DELTAX/6/,DELTAY/6/ 

BLACK=0 

BLUB=1 

RED=2 

MAGENT=3 

GREEN=4 

CYAN=5 

YELLOW-6 

WHITE=7 

COLOR=GKEEN 

IXSCLE=3 

ESC=27 

WRITE ( 5,499) ESC, ESC 
499  FORMAK '+' ,A1, 'PrTMl' ,A1, '®' ) 

C  !SET   TERMINAL   TO   ANSII    MODE 

TYPE   500 

500  FORMAT( IX, 'INPUT  FILE  NAME  ') 
ACCEPT  501,(FNAM(I) ,1-1 ,15) 

501  F0RMAT(15A1) 

OPEN (UNIT=1,NAME=FNAM,F0RM=' UNFORMATTED'  ,TYPE='OLD'  ) 
READ( 1) ,(DES(I) ,1=1 ,40) ,NRUN,NCNTS 
+,( (IY( J,K) ,K=1,NCNTS) ,J=1,NRUN) 
CL0SE(UNIT=1) 
TYPE   5 

5  F0RMAT( '       REVERSE   THE    SLICES?    (1-YES,    0-NO)     ') 
ACCEPT   6,NSLICE 

6  FORMAT (15) 
WRITE(5,61) 

61  PORMAT('  SHIFT  SUCCESSIVE  SLICES  (1  =  LEFT,  -1  =  RIGHT)  ' $) 
READ  (5,62)ISHFT 

62  FORMAT (13) 
C 

c 

WRITE(5,621) 
621  FOBMATC    WHAT    ARE  THE   VALUES    FOR  DELTAX,    DELTAY,    IXSCLE    [I]     '$) 
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READ(5,622)DELTAX,DEI,TAY,IXSCLE 

622  FORMAT(3I4) 
C 

C 

WRITE(5,623) 

623  FORMAT ('       WHAT    IS   CROSS    HATCH   INTERVAL    [I]     '$) 
READ(5,624)NXHTC 

624  FORMAT (13) 
C 

C 

DELTAX=DELTAX* ISHFT 
171      CALL  GGON 

CALL  GGINIT 

CALL  GGAXIS(0, 0,767, 479) 

CALL  GGERA 
C        CALL  GGBOX(7,0,0,767,479) 
C 

C        SET  UP  DATA  POINTS  FOR  X  AXIS 
C 

DO  19  NP0INT=1,NCNTS 

IX ( NPOINT ) =NPOINT* IXSCLE 
19       CONTINUE 

nline=1 

DO  11  NP0INT=1,NCNTS 

I0UTY( NPOINT )=IY( lABS ( NLINE-NSLICE*NRUN) ,NP0INT)/4 
1 1      CONTINUE 

DO  20  NLINE=1,NRUN,1 

DO  10  NP0INT=1,NCNTS 

I0UT(NP0INT)=IY( IABS(NLINE-NSLICE*NRUN) ,NP0INT)/4 
10      CONTINUE 

CALL  GG3DX( COLOR, IX, IOUT,IOUTY,NCNTS,NLINE,DELTAX,DELTAY,NXHTC) 

DO  30  NP0INT=1 ,NCNTS 

I0UTY( NPOINT ) =IOUT ( NPOINT ) 
30       CONTINUE 
20      CONTINUE 

CALL  AXIS ( COLOR, NRUN,NCNTS,DELTAX,DELTAY, IXSCLE) 
C 

c 

CALL  GGOFF 

WRITE( 5,2000) ESC, FN AM 

2000  format; '+' ,A1, ' [H'/'  ',20A1) 
WRITE(5,2001)ESC 

2001  F0RMAT( '+' ,A1 , ' [H  0-QUIT,  1 -SCREENDUMP ,  2-SCRDMP  NO  LABEL  '$) 
READ(5,2002)IANS 

2002  F0RMAT(I2) 
IF(IANS.EQ.0)GOTO2100 
WRITE(5,2004)ESC 

2004  FORMAT( '+■ ,A1, ' [H' ,80X) 
IP  (IANS.NE.2)G0T0  20035 
WRITE(5,2005)ESC 

2005  FORMAT( '+'A1, ' [H' ,80X/80X) 
20035    CALL  GGON 

WRITE(5,2003) 


215 


2003     FORMATC   S(H)') 

CALL  GGOFF 
2100     END 
C 
C 

c 
c 

SUBROUTINE  GG3DX( COLOR, IX, lY, lYX, NPNTS ,N,DELTAX,DELTAY,NXHTC) 

BYTE  COLOR 

DIMENSION  IX( 1) ,IY( 1) ,MASK( 1) ,IYX( 1) 

COMMON  /lAREA/MASK 

INTEGER  DELTAY,DELTAX 

IXOFP=200 

IF(DELTAX.LT.0)IXOFF=50 

IF(N.NE.1)GOT010 
C 
C        SET  UP  MASK  FOR  FIRST  SLICE 

DO  5  1=1 ,NPNTS 

MASK( IX( I)+( IXOFF-N*DELTAX) )=IY( I)+N*DELTAY 

5        CONTINUE 

C 

C 

1 0       CONTINUE 

DO  20  1=1 ,NPNTS 
IXOrjT=IX(  I)+(  IXOFF-N*DELT-AX) 
IYOUT=IY(  I)+N*DELTAY+20 
IF(IYOUT.GE.MASK( IXOUT) ) GOTO50 
GOTO  20 
5  0       MASK ( IXOUT ) =IYOUT 

CALL  GGPLT(COLOR, IXOUT, IY0UT,1) 
IF(N.LE.1)GOTO20 
IF(I.EQ.1)G0T019 
IF(IMOD(I,NXHTC).NE.0)GOTO20 

19  IX2=IX0UT+DELTAX 
IY2=IYX( I)+(N-1)*DELTAY+20 
IF(DELTAX.LT.O)GOT0190 
IF(IY2.LT.MASK(IX2))GOTO20 

190  CALL  GGVEC(COLOR,IX2,IY2) 

20  CONTINUE 
GOTO  33 

C        IF(N.EQ.1)GOT033 

C        M1=MASK(IX0FF-(N-1)*DELTAX) 

C        DO  33  1=1  ,DELTAX 

C        MASK(IX0FP-(N*DELTAX)+I)=M1 

33       CONTINUE 

RETURN 

END 


SUBROUTINE  AXIS ( COLOR, NLINE,NPNTS ,DELTAX ,DELTAY ,IXSCLE) 

BYTE  COLOR 

INTEGER  DELTAX,DELTAY,MASK( 1) 

INTEGER  X0,Y0,X1,Y1,X2,Y2,X3,Y3,X4,Y4,XS,Y5,X6,Y6,XORG,YORG 

INTEGER  X7,Y7 
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COMMON   /I AREA/MASK 
IXOFF=200 

IF(DELTAX.LT.O) IXOPP=50 
ISIGN=DELTAX/IABS( DELTAX ) 


X0=IXOFF 
Y0=20 


X1=IX0FF+IXSCLE*NPNTS 
Y1=20 


X2=IX0FF-NLINE*DELTAX 
Y2=2 0+NLINE*DELTAY 


X3=X2 

Y3=Y2+250 
C 

X4=X2+NPNTS*IXSCLE 

Y4=Y2 
C 

X5=X0 

Y5=Y0+250 
C 

X6=X4 

Y6=Y3 
C 

X0RG=X2 

Y0RG=Y2 
C 

X7=X1 

Y7=Y5 
C 

CALL   GGPLT(C0L0R,X6,Y6,1) 

IF ( DELTAX . LT . 0 ) GOTO 1 5 

IYTEMP=MASK(X4) 

CALL   GGVEC(C0L0R,X6,IYTEMP) 

GOTO 16 

15  CALL  GGVEC(C0L0R,X4,Y4) 
CALL   GGVEC(C0L0R,X1,Y1) 

16  CALL   GGPLT( COLOR, XI, Y1,1) 
CALL   GGVEC(COLOR,X0,Y0) 

IF    (DELTAX.GT.O)CALL   GGVEC(C0LOR,X2 , Y2 ) 

CALL   GGPLT(COLOR,X2,Y2, 1) 
161  IF(DELTAX.GT.O)GOT0191 

IYTEMP=MASK(X2) 

CALL   GGPLT(C0L0R,X2,IYTEMP,1) 
191  CALL   GGVEC(COLOR,X3,Y3) 

IF(ISIGN.GT.0)GOTO200 
C        SURFACE  FOR  LEFT  SHIFT 

CALL  GGPLT(COLOR,X3,Y3,1) 

CALL   GGVEC( COLOR, X5,y5) 

CALL  GGVEC(COLOR,X0,Y0) 
199  GOTO201 
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C  SURFACE    FOR   RIGHT    SHIFT 

200  CALL   GGPLT(COLOR,X6,Y6,1) 
CALL  GGVEC( COLOR, X7,Y7) 
CALL  GGVEC(C0L0R,X1 ,Y1) 

201  CONTINUE 
C 

C        VERTICAL  AXIS  TICS 

C 

DO  20  1=0,10 

IY0={I*2S)+V2 

IF ( DELTAX . LT . 0 ) GOTO 1 9 

CALL  GGPLT(COLOR,X2,IY0,1) 

CALL  GGVEC(COLOR,X2-6,IY0) 

GOT020 

19  CALL  GGPLT(COLOR,X4,IY0,1) 
CALL  GGVEC(COLOR,X4+6,IY0) 

20  CONTINUE 
C 

C        HORIZONTAL  AXIS  TICS 
C 

DO  25  1=0,10 

IX0=I*NPNTS*IXSCLE/10+X0 

CALL  GGPLT(COLOR,IX0,Y0, 1) 

CALL  GGVEC(COLOR,IX0,YQ-6) 
25       CONTINUE 
C 

C        ANGLE  AXIS  TICS 
C 

ZLINE=NLINE 

DO  35  ZI=0.,ZLINE,ZLINE/10. 

I F ( DELTAX . LT . 0 ) GOT027 

IXO=XO-DELTAX*ZI 

IYO=YO+DELTAY*ZI 

GOT023 

27  IX0=X1-DELTAX*ZI 
IY0=Y1+DELTAY*ZI 

28  CONTINUE 
IC4=ISIGN*6 

CALL  GGPLT(COLOR,IX0,IY0,1) 
CALL  GGVEC( COLOR, IX0-IC4,IY0-6) 

35       CONTINUE 

C 

C        BACK  AXIS  LINE 

C 

C        CALL  GGPLT( COLOR, X2,Y2,1) 

C        CALL  GGVEC( COLOR, X4,Y4) 

C 

C        TOP  AXIS  LINE 

C 

CALL  GGPLT( COLOR, X3,Y3,1) 
CALL  GGVEC( COLOR, X6,Y6) 
RETURN 
END 
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1  '  PROGRAM  MEASURES  WRITTEN  BY  JOHN  R.  RICHARDSON 

2  VERSI0N=2! 

3  '  VERSION  1.0  BASELINE  SET  6/1/85 

5  '  MAIN  PROGRAM  BEGINS  AT  LINE  1000 

6  '  VERSION  1.5  6/3/85 

7  '  CLEANED  UP  OLD  FORTRAN  CODE,  ADDED  DOUBLE  OUTPUT  FILE  MODE 

8  •  SAVE  TRUE  DIGITIZER  VALUES  FOR  PLOTTER  FILE  OUTPUT 

9  '  CALCULATE  AREAS  BASED  ON  SCALED  DATA 

10  '  VERSION  1.6  6/7/35 

11  '  FIXED  ERROR  IP  NO  FILES  OF  B  DRIVE, 

12  '  CHANGED  DATA  ARRANGEMENT  IN  OUTPUT  FILES 

13  '  VERSION  1.7  6/24/85 

14  '  ADDED  ERROR  OUTPUT  ROUTINE  FOR  ERRORS  OTHER  THAN  NO  FILES 

15  '  VERSION  2.0  ADDED  SCALE3  SUB  FOR  DIFFERENT  XSCALE  AND  YSCALE 

16  '  7/2/85 

99  '*****♦♦***««♦*«»♦.«««*♦».»»*«**««♦,,,*»»«♦«,«», »»»,,,»4,»,,,,»,,, 

100  ■  SUBROUTINE  DIGINI  (LINE  3000-3490)  OPENS  DIGITIZER 
110  '   AND  SETS  INITIAL   PARAMETERS  FOR  PROGRAM 

130  '  SUBROUTINE  STREAM  MODE  (3800-3899)  TURNS  ON  STREAM  MODE 
ISO  '  SUBROUTINE  POINT  MODE  (3900-3999)  TURNS  ON  POINT  MODE 
170  '  SUBROUTINE  FILE  HANDLER  (4000-4220)  OPENS  DATA  FILE  FOR  OUTPUT 
190  '  SUBROUTINE  SCALE2  (5000-5610)  HANDLES  SETTING  UP  USER 
200  '     COORDINATES  AND  ORIGINS  AND  OFFSETS 

220  '  SUBROUTINE  DELAY  (6000-6010)  ARE  TIMING  ROUTINES  THAT 
230  '   MAY  BE  NEEDED  FOR   SENDING  SETUP  INFORMATION  TO  DIGITIZER 
250  '  SUB  INPUT   (8000-8070)  GETS  DATA  SENT  FROM  DIGITIZER 
270  '  SUB  DIGURU  (9000-9070)  SCALES  DIGITIZER  INPUT  TO  REAL  WORLD 
290  •     COORDINATES 

300  '  SUBROUTINE  AREAP  (1660-2030)  GETS  INPUT  POINTS  FOR  AN  AREA 
330  '  SUBROUTINE  AREAX  (2070-2580)  CALCULATES  THE  AREA 
350  '  SUBROUTINE  PERIX  (2610-2810)  CALCULATE  THE  CLOSED  AND  OPEN 
370  '     PERIMETERS  FROM  A  SET  OF  GIVEN  POINTS 

1000  '*••*♦************♦*♦****♦******♦*♦***«*♦♦*♦**«♦•*****»♦♦**«»** 
1010  '*************************  MAIN  PROGRAM  START  ***************** 
1020  '*****************************«*************••*»«♦»**»«***»*«*♦ 
1040  ON  ERROR  GOTO  20000 
1060  CLS: PRINT: PRINT: PRINT: PRINT 

1070  PRINT  "AREA  MEASUREMENT  PROGRAM  VERSION  ", -VERSION 
1120  GOSUB  4000:'  CALL  FILE  HANDLER 
1130  GOSUB  3000:'  CALL  DIGINI 
1140  GOSUB  5000:'  CALL  SCALE2 
1260  'CONTINUE 

1270  GOSUB  1660:'  CALL  AREAP  ( X,  Y,  AREA,NPOINT,IERR,RESOL) 
1280   IF  (IERR=0)  THEN  GOTO  1290  ELSE  PRINT  "ERROR  ";IERR;" 
HAS  OCCURRED  NOT  ENOUGH  POINTS  FOR  AN  AREA":GOTO  1270 
1290  ' 
1340  GOSUB  2070:'  CALL  AREAX ( X,Y, AREA, NPOINT) 

1350  GOSUB  2610:'  CALL  PERIX( X ,Y,PERI1 ,PERI2 , NPOINT) 

1351  AREA=AREAIO 


1370  PRINT    "    THE   MEASURED   AREA    IS    [" 

1380  PRINT    "      CLOSED   PERIMETER= 

1390  PRINT    "      OPEN        PERIMETER=  " 

1400  PER10UT=PERI2 

1410  PRINT   BELLS  i"    KEEP   THIS   AREA  OR  RE-MEASURE    [K   or    R] 


AREAIO;"]     FOR   ";NP0INT;"    POINTS" 

PERI  2 

PERIl 
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1420  INPUT  ANS$ 

1440   IF  (ANS$  ="R")'30TO  1260 

1480  PRINT   #2,XT(1),YT(1),PENDUM;AREA;PERI0UT;NPOINT 

1481  PRINT    #3,AREA;PERI0UT 
1490    PRINT   #2,XT( 1),yT( 1),PENUP 
1500    FOR   1=2   TO   NPOINT 

1510    PRINT   #2,XT(I),YT(I),PENDWN 

1520  NEXT  I: '200       CONTINUE 

1525  IF  WET5=99  THEN  1590 

1530  PRINT  #2,XT( 1) ,yT( 1) ,PENDWN 

1590  PRINT  "  DO  YOU  WANT  TO  INPUT  ANOTHER  AREA  (Y  OR  N)  "  ,- 

1600  INPUT  ANS$ 

1620   IF  (ANS$  =  "Y")  GOTO  1260 

1630   CLOSE-PRINT  "TYPE  SYSTEM  TO  EXIT  FROM  BASIC  OR  RUN  TO  RERUN": END 

1640  .*..**».*•♦•♦*♦••**«*♦*•*********♦**********•******************* 

1641  i»*»«**»***«****«******«**    END   OF  MAIN   PROGRAM    ***************** 

1642  .*..♦*..♦•..*♦♦»•♦****•♦*•*♦************************************ 
1660    PRINT    "SUBROUTINE   AREAP" : IERR=0 

:'    SUBROUTINE    AREAP(X,Y,  AREA, NPOINT  ,IERR,RESOL) 

1700    IERR=0 

1710  NPOINT=0 

1720   AREA=0 

1740  PRINT  "  ENTER  FIRST  POINT  BY  PRESSING  THE  'V  KEY   " 

1745  PRINT  "  ENTER  REMAINING  POINTS  BY  PRESSING  ANY  KEY  BUT  '2'" 

1750  PRINT  "  THEN  QUIT  ENTERING  POINTS  BY  PRESSING  '2'  " 

1760  GOSUB  9000:'  CALL  DIGURU  (XIN, YIN , CODE) 

1770   XOLD=XIN 

1780   YOLD=YIN 

1790   IF  (CODES  <>  "1")  GOTO  1760 

1800  NP0INT=1 

1810  X( NPOINT )=X0LD 

1811  XT(NPOINT)=XTRUE 

1820  Y( NPOINT )=YOLD 

1821  YT(NPOINT)=YTRUE 

1830    IF  (CODES  =  "2")  GOTO  1950 

1840  GOSUB  9000: '200       CALL  DIGURU  ( XIN, YIN, CODE) 

1850    IF  (CODES  =  "2"  )  GOTO  1950 

1870    NP0INT=NP0INT+1 

1880  X(NPOINT)=XIN 

1881  XT ( NPOINT )=XTRUE 

1890  Y( NPOINT )= YIN 

1891  YT( NPOINT) =YTRUE 
1900    XOLD=XIN 

1910    YOLD=YIN 

1930  PRINT  XIN, YIN, "CODE  ="; CODES : BEEP 

1940    GOTO  1840 

1950  IP  (NPOINT  <  3)  THEN  lERR  =1:  '300 

1970    IF  (lERR  <>   0)THEN  RETURN 

2020    IERR!=0 

2025  'GOSUB  3900 

2030    RETURN 

2070  PRINT  "  SUBROUTINE  AREAX( X,Y,AREAIO, NPOINT)" 

2130  PRINT    "  ****  DIGITIZER  AREA  CALCULATION  *♦*♦" 

2140  NUMPNT=0: '200 


220 


2150      A1=0! 

2160      A2=0! 

2170      AREAIOO! 

2180      NUMPNT=NUMPNT+ 1 

2190    'C    READ    FIRST    PAIR 

2200         XF=X(NUMPNT) 

2210         !fF=y(NUMPNT) 

2220  XP=XF 

2230  YP=YF 

2240    ' 300  CONTINUE 

2250    'C 

2260   NUMPNT=NUMPNT+1 

2270  XC=X(NUMPNT) 

2280    YC=Y(NUMPNT) 

2290        A1=A1+(XP*YC) 

2300         XP=XC 

2310    IF(NUMPNT  =   NPOINT)GOTO  2330 

2320   GOTO  2240 

2330    IF(NUMPNT  >   2)G0T0  2370:'400 

2340  '       WRITE(LtJNO,401) 

2350  PRINT   "  7N0T  ENOUGH  DATA  POINTS  FOR  AREA  CALCULATION" 

2360         RETURN 

2370    A1=AH-(XC*YF)  : '450 

2380    NUMPNT=1 

2390    XF=X(NUMPNT) 

2400    YP=Y(NUMPNT) 

2410         XP=XF 

2420         YP=YF 

2430  NUMPNT=NUMPNT+1: '500 

2440    yC=Y(NUMPNT) 

2450    XC=X(NUMPNT) 

2460         A2=A2+(YP*XC) 

2470         YP=YC 

2480    IP  (NUMPNT  =   NPOINT)GOT0  2500 

2490         GOTO  2430 

2500   A2=A2+( YC*XP) :'600 

2510         AREAI0=ABS((A1-A2)*.5) 

2520         RETURN 

2530  '700    AREAIO=0. 

2540  '       RETURN 

2550  '800        CONTINUE 

2560  '       AREAlO=0. 

2570         RETURN 

2580  '       END 

2590  'C 

2600  'C 

2610  PRINT  "'  SUBROUTINE  PERIX(X, Y,PERI 1 ,PERI2 ,NPOINT) " 

2620  '   DIMENSION  X(1),y(1) 

2630    PER1 1=0! 

2640    PERI2=0! 

2650    NUMPNT=1 

2660    XIN=X(1) 

2670    YIN=y(1) 

2680    XL=XIN 
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2690 
2700 
2710 
2720 
2730 
2740 
2750 
2760 
2770 
2780 
2790 
2800 
2810 
2820 
2830 
3000 
3001 
3002 
3003 

3005 
3006 
3007 
3008 
3009 
3010 
3011 
3020 
3025 
3030 
3031 
3040 
3041 
3050 
3051 
3060 
3061 
3070 
3071 
3080 
3031 
3090 
3091 
3100 
3110 
3111 
3310 
3320 
3330 
3340 
3350 
3360 
3370 
3380 


YL=YIN 

troMPNT=NUMPNT+1 : ' 100 

XN=X(NUMPNT) 

YN=y(llUha>NT) 

PERI1=PER11+FNRDIST(XN,YN,XL,YL) 

XL=XN 

YL=™ 

IF  (NUMPNT  =   NP0INT)GOT0  2780 

GOTO  2700 

PERI2=FNRDIST    (XIN,yiN,XL, YL) : ' 300 
PER12=PER12+PER1 1 
RETURN 
END 


PRINT    "SUBROUTINE   DIGINI"'  SUBROUTINE   DIGINI 

LCB=13 

DEF  FNRDIST(X1 , Y1 ,X2 ,y2)=ABS( { X1-X2 )©2+( Y1-Y2 )®2)©. 5 

DEF  FNANGLER(X1,Y1,X2,Y2)=ATN( ( Y2-Y1 )/(X2-X1 ) ) 

+(SGN(ABS(X2-X1) ) -SGN( ( X2-X1 ) ) )* 1.570796 

PRINT  "INITIALIZATION  SEQUENCE  FOR  DIGITIZER" : BEEP 

PRINT"  PLEASE  MAKE  SURE  DIGITIZER  IS  ON" 

INPUT  "  HIT  RETURN  MHEN  READY  ",DUM$ 

DIM  X{ 1000) ,Y( 1000) ,XT( 1000) ,YT( 1000) 

PENUP=3 : PENDWN=2 : PENDUM=6 : BELL$=CHR$ ( 7 ) 

OPEN"COM1:9600,E,7,2,RS,CS,DS,CD"  AS  # 1 : ' OPEN  AUX  PORT  FOR  I/O 

CLS 

LD$="#]":QT$="/" 

PRINT  #1,"#]L": PRINT  "  DIGITIZER  BEING  RESET  ":GOSUB  6005 

PRINT  #1,"#] {":'SET  RESOLUTION  TO  .001 

GOSUB  6001 

PRINT  #1,"#] 6": 'SET  RATE  TO  2/SEC 

GOSUB  6001 

02": 'SET  INCREMENT  TO  .01 


PRINT  #1,"#] 

GOSUB  6001 

'PRINT  #1,"# 

'GOSUB  6001 

'PRINT  #1,"#]  J": 

'GOSUB  6001 

PRINT  #1,"#19":' 

GOSUB  6001 

PRINT  #1,"#] >": ' 

GOSUB  6001 

PRINT  #1,"#] I": ' 

PRINT  #1, "#]/":' 

GOSUB  6001 

XSCALE= 1 ! 
XOPF=0! :• 
YSCALE=1! 
YOFF=0! : 
ANGLE=0! : 

'       XROUND=0. 

'       YROUND=0. 
UXSCAL= 1 1 


M" : '  TURN  ON  INCREMENTAL  MODE 


TURN  ON  STREAM  MODE 

SET  SERIAL  TAG  AS  LAST  CHARACTER 

SET  NO  FIELD  DELIMITERS 

RESET  FOR  POINT  MODE  FOR  BEGINNING  SETUP 
SEND  END  OF  REMOTE  FORMATTING 

I ' USER  X  AXIS  SCALE  FACTOR 

USER  X  AXIS  OFFSET 
I ' USER  Y  AXIS  SCALE  FACTOR 

'USER  Y  AXIS  OFFSET 
'  USER  SKEW  CORRECTION  FACTOR 


'USER  PLOTTER  X  SCALE  FACTOR 
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3390         aYSCAL=  1 !  :  '  USER  PLOTTER  Y  SCALE  FACTOR 

3400         UROT=0!   :  '  USER  PLOTTER  ROTATION  ANGLE 

3410  '       ARRLEN=.15 

3420  ■       ARRWID=.07 

3430  '       ARROFP=.03 

3440  '       IARRTy=3 

3450  '       PDLEN=.1 

3460  '       PULEN=.05 

3470         XLAST=0  !  :  '  LAST  X  COORD  CALCULATED 

3480         yLAST=0! :'LAST  Y  COORD  CALCULATED 

3490  RETURN 

3500  'CLOSE  FILE  THEN  REOPEN  IT 

3510  CLOSE  #1 

3520  CPEN"C0M1 :9600,E,7,2,RS,CS,DS,CD"  AS  #1:'0PEN  AUX  PORT  FOR  I/O 

3530  RETURN 

3800  PRINT  "STREAM  MODE  SUBROUTINE" 

3810  PRINT  #1,"#]J":'  TURN  ON  STREAM  MODE 

3811  GOSUB  6001 

3820  PRINT  #1,"#1M":'  TURN  ON  INCREMENTAL  MODE 

3821  GOSUB  6001 

3899  RETURN 

3900  PRINT  "  SUBROUTINE  FOR  SETTING  POINT  MODE" 

3910  PRINT  #1,"#]I":'  RESET  FOR  POINT  MODE  FOR  BEGINNING  SETUP 

3911  GOSUB  6001 

3999  RETURN 

4000  PRINT  "FILE  HANDLING  SUBROUTINE":"  SUBROUTINE  FOR  FILE  HANDLING 

4004  PRINT  "CURRENT  DATA  FILES:  ":PRINT 

4005  FILES  "B:*.*" 

4006  PRINT: PRINT: PRINT: 
4010  'PRINT  "FILE  HANDLING" 

4020  PRINT  "WOULD  YOU  LIKE  TO  OPEN  A  NEW  FILE  OR  APPEND  TO  AN  EXISTING" 
4030  INPUT  "    FILE  ( 1-NEW,  2-OLD) " ,FILEMODE 
4040  IF  FILEMODE  <1  OR  FILEMODE  >2  THEN  4020 
4050  IF  FILEMODE  =2  THEN  4100 

4060  INPUT  "WHAT  IS  THE  NAME  FOR  THE  FILE  (1-8  CHARACTERS)  "[FILENAMES 

4061  IF  LEN  (FILENAME$)>8  THEN  4060 

4062  IF  INSTR( FILENAMES ,":" )  <>0  THEN  PRINT  "INPUT  FILENAME 
ONLY  WITHOUT  DRIVE  SPECIFIER" :GOTO  4000 

4070  NTEMP=INSTR( FILENAMES,".") 

4075  IF  NTEMP=0  THEN  4085 

4078  NLEN=LEN( FILENAMES) 

4080  FILENAMES=MIDS ( FILENAMES , 1 , NTEMP- 1 ) 

4085  FILENAMES="B:"+FILENAMES 

4090  OPEN  "0", 2, FILENAMES*". DAT" 

4091  OPEN  "0",3,FILENAMES+".PPN" 

4093  INPUT  "WHAT  IS  DESCRIPTION  OF  THIS  DATA  SET" ,DESCS 

4094  PRINT  #3,CHRS(34)+DESCS 

4095  GOTO  4220 

4100  'OPEN  FOR  APPEND 

4110  INPUT  "WHAT  IS  THE  NAME  OF  THE  EXISTING  FILE  ", FILENAMES 

4120  IF  LEN( FILENAMES )>8  THEN  4110 

4130  IF  INSTR(FILENAMES,":")  <>0  THEN  PRINT  "INPUT  FILENAME 

ONLY  WITHOUT  DRIVE  SPECIFIER" :GOTO  4100 
4140  NTEMP=INSTR( FILENAMES,".") 
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4150  IF  NTEMP=0  THEN  4200 

4160  NLEN=LEN( FILENAMES) 

4170  FILENAME$=MID$ ( FILENAMES , 1 ,NTEMP- 1 ) 

4200  FILENAME$="B:"+FILENAME5 

4202  OPEN  "I",3,FILENAME$+".PRN" 

4204  LINE  INPUT  #3 ,DUM$ :PRINT: PRINT  DUMS : PRINT : PRINT 

4206  CLOSE  #3 

4210  0PEN"A",2,FILENAME$+".DAT" 

4211  0PEN"A",3,FILENAMES+".PRN" 
4220  'INPUT  "BASIN  NUMBER  =  ",HET1 
4250  RETURN 

5000  'PRINT'-SUBROUTINE  SCALE3"  :  '       SUBROUTINE  SCALES 

5030  GOSUB  3900:  •  GO  SET  POINT  MODE  FIRST  1 

5040  PRINT  "  •*•♦**  DIGITIZER  THREE-POINT  SCALING  •*•***" 

5050  PRINT  :PRINT:PRINT  :BEEP:GOSUB  6001 

5060  PRINT  "   DIGITIZE  THE  ORIGIN  OF  THE  GRAPH  >>":BEEP 

5070  GOSUB  8000:XORG=XIN!YORG=YIN:'       CALL  DIGDRn(XORG, YORG,IBTN) 

5080     IF(VAL(CODE$)=10)GOTO  5530 

5100  PRINT  "  DIGITIZE  ANY  OTHER  KNOWN" 

5110  PRINT  "   POINT  ON  THE  SAME  HORIZONTAL  (X-AXIS)  LINE>>" 

5120  GOSUB  8000:XHZ=XIN:YDUMM=YIN: '       CALL  DIGDRU(XHZ, YHZ ,IBTN) 

5130    IP(VAL(CODE$)=   lOjTHEN  RETURN 

5140         XSCALB=1! 

5150         XOFF=0! 

5160         YSCALE=1! 

5170         YOFF=0! 

5180         ANGLE=0! 

5190         XROUND=0! 

5200         YROaND=0! 

5210         XD=FNRDIST(XORG,YORG,XHZ,YDUMM) 

5220         IF(XD  =   0!)THEN  RETURN :' GOSUB  3800 :  RETURN: ' STREAM  MODE 

5240  'PRINT  "  TWO-POINT  X  DISTANCE  IN  DIGITIZER  REAL  UNITS:  ";XD 

5242  PRINT  "DIGITIZE  A  THIRD  KNOWN  POINT  ON  THE  VERTICAL  (Y-AXIS)  LINE" 

5244  GOSUB  8000:XDUMM=XIN:YHZ=YIN 

5246  ANG1=FNANGLER(X0RG,Y0RG,XHZ,YDUMM) 

5247  RNG2=FNANGLER(XORG,Y0RG,XDUMM,YHZ) 

5248  YANG=1.570796-(ANG2-ANG1) 

5249  YD=PNRDIST( X0RG,Y0RG,XDUMM, YHZ) *COS( YANG) 

5250  IP  YD=0  THEN  RETURN 

5252  PRINT  "X-DISTANCE=  ";XD;"   Y-DISTANCE=  ";YD 

5260  PRINT  "  ENTER  USER  COORDINATES" 

5270  PRINT  "  OF  THE  FIRST  POINT  (REAL)  [0.0,0.01:  " ; 

5280  INPUT  X1U,y1U 

5290  '5      FORMAT(2F10.0) 

5300         XDEF=XD+X1U:YDEF=YD+Y1U 

5310  •       WRITE(LUN0,6)XDEF 

5320  PRINT   "  ENTER  USER  X  COORDINATE" 

5330  PRINT  "  OP  THE  SECOND  POINT  (REAL)  " ;XDEF 

5340  INPUT  X2U 

5350  IP(X2U  =   0! )THEN  X2U=XDEF 

5351  PRINT  "  ENTER  USER  Y  COORDINATE" 

5352  PRINT  "  OF  THE  THIRD  POINT  (REAL)  ";YDEP 

5353  INPUT  Y3U 

5360         XU=X2U-X1U 
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5365         YU=Y3U-ylU 

5370         ANGLE=ANG1 

5380         IF(XU  <>  0!)THEN  XSCALE=XU/XD 

5390         IF(TO  <>  0)THEN  YSCALE=yU/YD 

5400         IP(X1U  =   0!  AND  Y1U  =   0!)GOTO  5510 

5402         X 1 0=X  1 U/XSCALE :  Y 1  U=  Y1 U/ YSCALE 

5410         ANGLU=FNANGLER(01 ,01 ,X1U,Y1U) 

5420         ANGLU=ANGLE+ANGLU 

5430         DISTU=FNRDIST(0! ,0! ,X1U,Y1U) 

5440  ■        DISTX=D1STU/XSCALE 

5450  '       DISTY=DISTU/YSCALS 

5460         XROT=DISTU*COS(ANGLU) 

5470         yROT=DISTU*SIN(ANGLU) 

5480         XOFP=XORG-XI!OT 

5490         YOFP=YORG-YROT 

5500         GOTO  5530 

5510         XOFP=XORG: ' 100 

5520         YOPF=YORG 

5530  '200    HRITE(LUN0,2ai) 

5540  PRINT  "  ENTER  X-AXIS  ROUNDOFF  (REAL)  [0.0]:  "  ; 

5550  INPUT  XROUND 

5560  '202    FORMAT(F6.0) 

5570    '  HRITE(LUNO,203) 

5580    PRINT    "    ENTER  Y-AXIS    ROUNDOFF    (REAL)     [0.0]:    " ; 

5590    INPUT    YROUND 

5600    RETURN :' GOSUB    3800:RETURN    :' RESET    STREAM   MODE   FIRST! 

5610    '  END 

6000  '  TIMER  LOOPS 

6001  BEEP  :FOR  IDUM=1  TO  375  :NEXT  IDUM:PRINT  TIMES : RETURN :' 1  SEC 

6002  BEEP  :FOR  IDUM=1  TO  750  :NEXT  IDUM:PRINT  TIMES : RETURN :' 2  SEC 

6003  BEEP  :FOR  IDUM=1  TO  1125:NEXT  IDUM:PRINT  TIMES : RETURN: ' 3  SEC 
6005  BEEP  :FOR  IDU»=1  TO  ia75:NEXT  IDUM:PRIN'r  TIMES :  RETURN :' 5  SEC 
6010  BEEP  :FOR  IDUM=1  TO  3750:NEXT  IDUM:PRINT  TIMES :RETURN: • 1QSEC 
8000  REM  ++++++  GET  INPUT  FROM  COM  BUFFER  ++++++++ 

8010  WHILE  L0C(1)  <  LCB 

8020  HEND 

8030  'IF  L0F(1)  <  24  THEN  BEEP:BEEP:  BEEP 

8040  DZ$=INPUTS(LCB,#1) 

8041  'PRINT  DZS 

8050    XS=LEFTS(DZ$,5)     :    Y$=MIDS  (  DZS  ,6  ,5  )     :    CODES=    rilDS  (  DZ$  ,  1 1  , 1 ) 

8060  XIN=VAL(XS)/1000    :    YIN=VAL( YS )/1 000 

8061  ' PRINT   XIN, YIN, CODES:BEEP 
8070    RETURN 

9000    '     SUBROUTINE   DIGURU( X, Y,IBTN) 

9010    GOSUB   8000 

9020    DISTU=P!'IRDIST(XOFF,yOFF,XIN,YIN) 

9030  ANGU=FNANGLER( XOFF , YOPP , XIN , YIN) -ANGLE 

9040  XIN=DISTU*COS(ANGU)*XSCALE 

9050  yiN=DISTU*SIN(ANGU)*YSCALE 

9055  XTRUE=XIN/XSCALE :  YTRUB=  YIN/XSCALE 

9070  RETURN 

20000  IP  ERR=53  AND  ERL=4005  THEN  PRINT  "NO  FILES  ON  B:":RESUME  4006 

20030  PRINT  "ERROR  NUMBER  ";ERR;"  HAS  OCCURRED  AT  LINE  ";ERL 
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PULSING  MODEL  GIGI  VERSION 

MODIFIED  TO  WRITE  DATA  PILE  FOR  3-D  GRAPHICS  PROGRAM 

PULSEGGM.FTN  VERSION 

BASELINE  MODEL  2.0  PREVIOUS  TO  AUG  29,1984 

VERS  2.1  CHANGED  FORMAT  OF  HARDCOPY  PRINTOUT  OF  VARIABLES 

VERS  2.11  ADDED  TOTALS  AND  DESCRIPTION  TO  OUTPUT  LIST 

VERS  2.12  ADDED  K8  TO  OUTPUT  LIST 

VERS  2.121  (1/24/85)  ADDED  TRACKING  COEFFICIENTS  ON  K2,K9,K11 
AND  CHANGED  FORMAT  11003  FOR  VARS(IVAL)  FROM  F6.1  TO  G15.4 
PROGRAM  PULSE 

BYTE  XTEXT(80) ,DES(40) ,YTEXT(80) 
DIMENSION  VARS(20),ALPHA(20) 
DIMENSION  FILE(3) 

REAL  M1,M2,M3,M4,M9,K10,X11,K12,K13 
REAL  K1,K2,K3,K4,K5,K6,K7,K8,K9,J,J0,JNORM 
EQUIVALENCE  ( VARS ( 1 ) ,K1 ) , 

VARS(2),K2), 

VARS(3) ,K3) , 

VARS(4) ,K4), 

VARS(5) ,KS) , 

VARS(6) ,K6), 

VARS (7) ,K7) , 

VARS(8) ,K8) , 

VARS(9) ,K9) , 

VARS( 10) ,K10), 

VARS(11),K11), 

VARS( 12) ,K12), 

VARS(13) ,K13), 

VARS( 14) ,XJ0INI), 

VARS( 15) ,21IC), 

VARS( 16) ,Q2IC), 

VARS( 17) ,Q3IC), 

VARS( 18) ,Q4IC) 
VIRTUAL  IY(6,25,180) 


DATA  FILE/3*'     '/ 

DATA  XTEXT/'  ' , 'P' , 'u ' , ' 1 ' , ' s ' , 'e ' , '  ' , 'M' 
'I','  ',67*0/ 
DATA  DES/40*0/ 
DATA  YTEXT/80*0/ 
DATA  ALPHA/' K1   ' , 'K2   ' , 'K3 

,'K8   ','K9   '.'KIO  •,'K11 

,'Q2IC' , 'QSIC' , 'Q4IC' ,2*' 
FT1(A,B)=ABS(AINT{A/B)-A/B) 

VERS=2.121  19/7/84    ;     1/24/85 

WRITE(5,100) 

FORMAT('       WHAT    IS   THE    DATA   FILE   FOR  THIS    MODEL    RUN    7') 
READ(5,101) (FILE(I) ,1=1,3) 
F0RMAT(3A4) 


K4      • 

,'K5      ' 

,'K6       ' 

'K7       • 

K12    • 

,'K13    ' 

, 'JOIN' 

'QIIC 

/ 
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MRITE(5,1011)PILE 

1011  FORMAT ( IX, 3A4) 
C                   WRITE(5,1016) 

C1016         FORMAT ( '       DO   YOU   WANT    HARDCOPY    (1-YES,    0-NO)     '$) 
C  READ(5,1017)ICOPy 

ICOPY=0 
C1017         F0RMAT(I2) 
C  WRITE    (5,1018) 

C1018         FORMATC       WHICH   Q   TO    SAVE    (  1 ,2  ,  3  ,4  ,5=JR,6=%POW  USED)'?) 
C  READ    (5,1019)I2SAV 

I2SAV= 1 
C1Q19         P0RMAT(I3) 
C  WRITE(5,1020) 

C1020         FORMAT (  '    DO  YOU   WANT   TO   PLOT   THE   GRAPHS    (1-YES,    0-NO)'$) 
C  READ(5,1021)IPLOT 

IPLOT=0 
C1021         FORMAT(II) 
C  WRITE(5,99) 

C99  FORMATC  HOW  LONG   TO    RUN?    ') 

C  READ(5,98)TIME 

C98  FORMAT (G6.0) 

CALL   ASSIGN( 1 ,FILE) 

READ( 1) El ,E2,E3,E4,E5,E6,E7 ,E8,E9,E10,E11 ,E12,B13,E14,E15, 
C  COEFFICIENTS    ********************************************** 

+NUM,K1,K2,K3,K4,KS,K6,K7,K8,K9,K10,K11 ,K12,K13 
+,XJ0INI,Q1 ,Q2,Q3,Q4 
C  INITIAL   CONDITIONS***************************************** 

CLOSE    (UNIT=1) 

Q1IC=Q1 

Q2IC=Q2 

Q3IC=Q3 

Q4IC=Q4 
(^******************************************************** 

C 
C 

WRITE(5,1012)K1  ,K10,K2,K11  ,.K3  ,K12  ,K4  ,K1  3  ,K5  ,XJOINI  ,K6  ,Q1  , 
+K7,Q2,K8,Q3,K9,Q4 

1012  FORMAT( 1X,'1-K1       ',G12.6,'  10-K10 ' ,G12 .6/ 
+1X,'2-K2       ',G12.6,'  11-K11' ,G12.6/ 
+1X,'3-K3       ',G12.6,'  12-K12' ,G12.6/ 
+1X,'4-K4       ',G12.6,'  13-K13' ,512.6/ 
+1X,'S-K5       ',G12.6,'  14-XJOINI' ,G12.6/ 
+1X,'6-K6       ',G12.6,'  15-Q1IC' ,G12.6/ 
+1X,'7-K7       ',G12.6,'  16-Q2IC' ,G12.6/ 
+1X,'8-K8       ',G12.6,'  17-Q3IC' ,G12.6/ 
+1X,'9-K9       ',G12.5,'  18-Q4IC' ,G12.6/ 

+1X,'  INPUT  VARIABLE   NUMBER  TO   VARY  =>    '$) 

READ(5,1013)IVAL 

1013  F0RMAT(I2) 

WRITE    (5,1014)ALPHA(IVAL)  ,VARS(IVAL) 

1014  FORMAT ( '       VARIABLE    ',A4,'    =    ',G12.6/ 
+'       HOW  MUCH   TO    INCREMENT?    '$) 

READ(5,1015)XINC 

1015  F0RMAT(G15.6) 
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C 
C 
C 

Q*4 ****************************************************** 

CALL  ASSIGN  (4,'TTO:') 

IF  (IPLOT.EQ.0)GOTO499 

CALL  GGON 

CALL  GGINIT 

CALL  GGAXIS(0, 0,767, 479) 

499  Q1IC=Q1 
Q2IC=Q2 
Q3IC=Q3 
Q4IC=Q4 
NSLICE=25 
NCNTS=150 
NRUN=0 

500  CONTINUE 

C  J0=XJ0INI+NRUN*4. 

IF(NRtlN.EQ.0)GOTO501 

VARS(rVAL)=VARS(IVAL)+XINC  [INCREMENT   VALUE  WE   ARE   VARYING 

501  NRUN=NRUN+1 
JO=XJOINI 

J=JO/( 1+K13*Q1*Q4)       IGIVE  JR  (J)  INITIAL  VALUE 
T=0. 
Q********************************************************** 

C        SET  K  VALUES  TO  TRACK  FOR  MULTIRUN  MODEL 


C 


K3=.1*K2 
K4=.9*K2 


K7=.1*K9 
K8=.9*K9 


K5=.1*K11 
K6=.9*K11 

Q* ********************************************************** 

Q1=Q1IC 

Q2=Q2IC 

Q3=Q3IC 

Q4=Q4IC 

EUSED=0.a 

M1=0. 

M2=0  . 

M3=0. 

M4=0  . 

R1=0. 

R2=0. 

R3=0. 

R4=0. 

R5=0. 

R6=0. 

R7=0. 

R8-0. 

R9=0. 
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R10=0. 

R11=0. 

R12=0. 

EUSED=0.0 

PAVAIL=0 . 0 
C 

C  HRITE(5,1081) 

C1081         FORMAT!'       SCALE   FACTOR  FOR  22 — 200.    OR    1000.  —  [  R]  ' ) 
C  READ(5,1082)SFACT 

C1082         FORMAT(G7.2) 
C  WRITE(5,108) 

C108  FORMAT('         WHAT    IS   THE   TIME   INTERVAL   DT    [  R]     ') 

C  READ(5,109)DT 

C109  FORMAT(G5.3) 

TIME=750. 

DT=.1 

SFACT=10D. 

NTIME>=TIME 

XDT=DT/10. 
C  START   OF   LOOP   *********************************************** 

IF    (IPLOT.EQ.0)GOTO5 

CALL   GGERA 

CALL   GGBOX(7,0,0,767,479) 

CALL  GGTEXT(7,626,475,XTEXT,1,0) 

CALL   GGBOX(7,620,452,767,479) 
C  PRINT    INITIAL  CONDITIONS************************************* 

5  T=T+DT 

C  RATE    EOUATIONS*********************************************** 

J=J0/(1+K13*Q1*Q4) 

POWUSE=100.*( J0-J)/J0 

PAVAIL=PAVAIL+JO 

EUSED=EUSED+JO-J 

R1=DT*K1*Q1*Q4*J 

R2=DT*K2*Q1 

R3=DT*K3*Q1 

R4=DT*K4*Q1 

R5=DT*K5*Q2 

R6=DT*K6*Q2 

R11=DT*K11*Q2 

R7=DT*K7*Q2*Q3*Q3 

R8=DT*K8*Q2*Q3*Q3 

R9=DT*K9*Q2*Q3*Q3 

R1 0=DT*K1 0*Q1*Q4* J 

R12=DT*K12*23 
C        LEVEL  EOUATIONS  ********************************************* 
1091    CONTINUE 

Q1=QH-R1-R2 

Q2=Q2+R3-R9-R11 

Q3=Q3+R5+R7-R12 

Q4=Q4+R4+R6+R1 2+R8-R1 0 

IF(FT1(T,1. ).GT.DT)GOT0    110 

IF(IPLOT.EQ.0)GOTO20000 

ITIME=T 

IXC=Q1/10. 
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CALL   GGPLT(4,IT1ME,IXC,1) 

IXC=Q2/SFACT 

CALL  GGPLT(3,ITIME,IXC,1) 

IXC=Q3/10.0 

CALL   GGPLT( 1,ITIME,IXC,1) 

IXC=Q4/10C. 

CALL  GGPLT(5,ITIME,IXC,1) 

IXC=(J0-J)*(250./J0) 

CALL  GGPLT(2,ITIME,IXC,1) 


! Q1    GREEN 

!  Q2    MAGENTA 

! Q3    BLOE 

!Q4    CYAN 

! POWER  USED    FED 


C 

C        FIND  WHICH  Q  TO  SAVE  IN  ARRAY 

20000    CONTINUE  !GOTO( 21000 ,22000 ,23000 ,24000 ,25000 ,26000) IQSAV 
C        GOTO110 

21000    IY( 1,NRUN,INT(T/5.)+1)=INT(Q1/2. ) 
C        GOT0 1 1 0 

22000    IY(2,NRUN,INT(T/5.)+1)=INT(Q2/20.) 
C        GOTO110 

23000    IY( 3,NRUN,INT(T/5. )+1 )=INT( Q3/2. ) 
C        GOT0110 

24000    IY(4,NRUN,INT(T/5. )+1 ) =INT( Q4/40 . ) 
C        GOTO  1 1 0 

25000    IY{5,NRUN,INT(T/S.)+1)=INT(( J0-J)*5. ) 
C        GOTO110 

26000    IY(6,NRUN,INT(T/5.)+1)=INT((POWUSE-80)*50) 
110     CONTINUE 

M1=AMAX1(M1,Q1) 

M2=AMAX1(M2,Q2) 

M3-AMAX1(M3,Q3) 

M4=AMAX1(M4,Q4) 

IF(T.LT.TIME)GOTO  5 
C 

ENCODE(80,11003,YTEXT)NRUN,ALPHA(IVAL)  ,VARS(IVAL) 
+,EnSED, 100*EUSED/PAVAIL 
11003    FORMAT (2X, 12, '  VARIABLE  ',A4,'  =  ',G15.4,'  POWER  USED  ', 
+G12.6,'  PPU:  ',G12.6) 

IF  (IPLOT.EQ.0)WRITE(4,11004)YTEXT 

FORMAT( 1X,80A1) 

IF  (IPLOT.EQ.DCALL  GGTEXT(  7  ,0  ,460  ,YTEXT,  1  ,0  ) 

IF  (ICOPY.EQ. 1)WRITE(3, 11001) 

FORMAT( ■+S(H)  '  ) 

IF(NRUN.LT.NSLICE)GOTO500 

IF  (IPLOT.EQ.0)GOTO24999 

CALL  GGERA 

CALL  GGOFF 

IQSAV=1 

CONTINUE 

CALL  ASSIGN ( 2 , ' PULSAV ' , 6 ) 

ENCODE( 40, 25001 ,DES) IQSAV, FILE 

F0RMAT(1X,'  TANK  Q',I1,'  FOR  DATA  FILE  ',3A4) 

WRITE(2)(DES(I) ,1=1 ,40) ,NRUN,NCNTS, 

( (IY( IQSAV, J1 ,K) ,K=1,NCNTS) ,J1=1,NRUN) 

CL0SE(UNIT=2) 

IQSAV=IQSAV+1 

IF  (IQSAV.LT.7)GOT024995 


11004 


11001 


24999 
24995 


25001 
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C 


WRITE(4,114) 
114  FORMAT(//'  Q1  Q2  Q3  Q4  TOTAL') 

WRITE(4,1121)Q1IC,Q2IC,Q3IC,Q41C,Q4IC+Q3IC+Q21C+Q1IC 
WRITE(4,1122)M1,M2,M3,M4 
WRITE(4,1123)Q1,Q2,Q3,Q4,Q4+Q3+Q2+Q1 

1121  FORMATC     INIT    '  ,4(  2X,G8  .2  )  ,2X,G12  .6) 

1122  FORMATC    MAX       '  ,4(  2X,G3.2)  ) 

1123  FORMATC    FINAL' ,4(  2X,G8.2  )  ,2X,G12  .6) 

WRITE(4,116)K1,K2,K3,K4,K5,K6,K7,K8,K9,K10,K11 ,K12 ,K13 , JO ,J 
116  F0RMAT(/1X, 'K1=    ',G12.6,'         K2=    ',G12.6,'  K3=    ',G12.6/ 

+'  K4=  ',G12.6,'  KS=  '.Gia.e,'  K6=  ',G12.6/ 
+'  K7=  ',G12.6,'  K8=  ',G12.6,'  K9=  ',012.6/ 
+■  K10=' ,G12.6,'  K11=' ,G12.6,'  K12=',G12.6/ 
+'    K13=' ,G12.6,'         J0=    ',G12.6,'         Jft=    ',G12.6) 

EUSED=EUSED*DT 

PAVAI]>PAVAIL*DT 

PPU= 1 0  0 . *EUSED/PAValL 

WRITE(4, 119) VERS, DT,EUSED,PAVAIL,PPU 

119  FORMAT(/'    PULSE    MODEL   VERS',F6.3,'    TIME   STEP(DT)    =    " ,F6.4/ 
+'    TOTAL   POWER   USED        =    ',G15.6,'    POWER   AVAILABLE  =',G15.6/ 
+'    PERCENT   POWER  USED   =    ',G15.6) 

WRITE(4,120)(FILE(I) ,1=1,3) ,E1 ,E2 ,E3 ,E4 ,E5 ,E6 ,E7 ,E8,E9,E10 , 
+E11,E12,E13,E14,E15 

120  FORMAT ( '    DATA   FILE  NAME   =    ' , 3A4 , IX, 15A4) 

I F ( ICOPY . EQ . 0 ) GOT0999 
CALL  GGON 
WRITE(3,11001) 
CALL  GGOFF 

END 


C 
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C        SURFACE  PULSING  MODEL  PROGRAM   3/24/83 

C        ADDITION  OF  CONSUMER  CEILING  TO  ALLOW  UP  TO  100  TOTAL 

C        CONSUMERS 

C 

C        DIFFUSION  ADDED  7/21/83  TO  NUTRIENT  TANK  Q4 

C 

C        VERS  3.01  ADDED  STARTING  CONDITION  TO  FILE  OUTPUT 

C        VERS  3.02  ADD  TIME  AND  DATE  TO  BEGINNING  OP  PROG 

C 

C 

PROGRAM  SURPUL 
C        Q1=PR0DUCER 
C        22=STORAGE  (PRODUCER) 
C        Q3=C0NSUMER 
C        Q4=NUTRIENTS 

DIMENSION  Q1( 12,12 ),Q4( 12,12) ,E(12, 12), Q3( 100) 

DIMENSION  Q4T( 12,12) 

DIMENSION  ETVPE(3) ,IX( 144) 

DIMENSION  Q2( 12,12) ,IXYZ( 100) 

INTEGER*4  ICNT(12,12) 

BYTE  TITLE( 10) ,ICON( 12,12) ,BUF1(9) ,BUF2(8) 

REAL  M,K1 ,K2,K3,K4,K5,K6,K7,K8,K9,K10,MTOT 

REAL  K11 ,K12,K13,J0,JR 

BYTE  ESC, TEXT(80) , COLOR, ICOLOR, CHAR 

INTEGER  X1(100),Y1(100) ,T1 ,T2 ,XTEMP,YTEMP 

FT1(A,3)=ABS(AINT(A/B)-A/3) 

IXY(I,J)=(I-1)*12+J 

DATA  TITLE/'D' , 'S' , 'P' , ' 1' , '0' , '0' , '  ','  ','  ','  '/ 

DATA   Q4T/144*0.0/  !  DK 

DATA   ICON/ 144*0/ 

DATA   BTYPE( 1)/'HIER'/ 

DATA   ETYPE(2)/'EVEN'/ 

DATA   ETYPE(3)/'RAND'/ 

ESC=27 

VERS=3.02 

CALL   TIME(BUF2) 

CALL   DATE(BUFI) 

WRITE( 5 ,5 ) ESC , ESC , ESC , ESC , TITLE, VERS ,BUF2 ,BUF1       13.0 

5  FORMAT( 1X,A1, 'PrTMl' ,A1, '»' ,A1, ' [2J' ,A1, ' tH' ,  !3.0 
&'    SURFACE   MODEL    ',10A1,'    — VERSION —    ',F5.2/                 ]3.0 
S1X,8A1,1X,9A1/                                                                                              13.02 
S'    DO   YOU   WANT   GRAPHICS    ON    (1-YES,    0-NO)     '$) 

READ    (5,6)IOFLAG 

6  FORMAT(II) 
C 

C 

WRITE    (5,7) 

7  FORMAT ( '    PLOTTING    INTERVAL   FOR  PRODUCER   AND   CONSUMER    [I]     '$) 
READ    (5,8)ITINT,ITINTC 

8  PORMAT( 13,13) 
WRITE( 5,808) 

808  FORMAT( '       HARDCOPY   AT    PLOTTING    INTERVAL    ( 1-YES, 0-NO) ' S) 

READ(5,8081)IPTR 
8081  F0RMAT(I3) 
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TINT=ITINT 
TINTC=ITINTC 
WRITE    (5,81) 
8082  CONTINUE 

81  FORHATC       HOW  LONG   TO    RUN?    [R]     '$) 
READ    (5,82)TTIME 

82  PORMAT(G6.0) 
WRITE(5,83) 

83  FORMAT ('    WHAT    IS    DT    [R]     '$) 
READ    (5,84)DT 

84  FORMAT(G10.6) 
XDT=DT/TINT 
XDTC=DT/TINTC 

WRITE(5,841)  IDK 

841  FORMAT( '  WHAT  IS  NUTRIENT  CONC.  OF  OUTER  NONREACTIVE ' /   !DK 
&  '    RING  (39000  IC;  0.0  TO  7  )  [R]  'S)  !  DK 

READ  (5,842)Q40IC  IDK 

842  FORMAT(G16.5)  IDK 
WRITE(5,85)  !DK 

35       FORMAT ( '   WHAT  IS  DIFFUSION  COEFFICIENT?  [  R]  '$)         ! DK 
READ  (5,36)DK  ! DK 

86       FORMAT (F8. 5)  IDK 

WRITE(5,9) 

9        FORMAT ('  INPUT  THE  SEARCH  LENGTH  FOR  PREDATOR  [I]  ',$) 
READ  ( 5 ,  1 1 ) N 

1 1  FORMAT ( 12 ) 
WRITE(5,91) 

91  F0RMAT(1X,'    FEEDING    AND   DOUBLING   THRESHOLD    [R,R]     'S) 
READ(5,92)PTHRSH,THRESH 

92  FORMAT (2G8. 2) 

WRITE(5,121)  13.0 

121  FORMAT ( IX, 'INPUT  ( 0-SUCCESSION;  1-STEADY  STATE)  [I]  '$)  !3.0 
READ(5,122)ISSUC  !3.0 

122  F0RMAT(I2)  !3.0 
WRITE(5,12) 

12  FORMAT( '  WHAT  ENERGY  TYPE  WOULD  YOU  LIKE'/ 
+   '1:  STD  INPUT'/ 

+   '2:  CONSTANT  INPUT'/ 

+   '  3:  RANDOM  INPUT' , 20X, ' ENERGY  TYPE  [I]  '$) 
READ(5,13)  lETYP 

13  PORMAT(II) 

C        IFdETYP.EQ.  1  )GOTO150 
WRITE(5,14) 

14  FORMAT('  WHAT  IS  THE  MEAN  VALUE  OF  ENERGY  '$) 
READ(5, 15)XMEAN 

1 5  FORMAT ( F5 . 2 ) 
150      CONTINUE 

GPP=0.  !GPP  COUNTER  (TOTAL) 

CNSUMP=0.        ICONSUMPTION  BY  CONSUMERS  (TOTAL) 

Q1TOT=0.         lAMOUNT  OF  PRODUCERS 

Q2TOT=0.         ! AMOUNT  OF  STORAGE 

Q3TOT=0.         [AMOUNT  OF  CONSUMERS 

Q4TOT=0.         ! AMOUNT  OF  NUTRIENTS 

PRDD=0.  I  TOTAL  PRODUCTION  (GPP+CNSUMP) 
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ETOT=0.  ! TOTAL   ENERGY    INPUT    (SUM  OF   MATRIX) 

TOTPOW=0.0  ! MEASURE   TOTAL   POWER  USED:    SUM  OF   EUSED 

T=0.  1 

T1=1  ! NUMBER  OF  CONSUMERS 

T2=T1 

Q1IC=1000.       !IC  CONDITION  FOR  Q1 

Q4IC=39000. 

Q2IC=1000. 

IF(ISSUC.EQ.0)GOTO147  !3.0 

e4IC=30000.      !IC  CONDITION  FOR  NUTRIENT  TANK  !3.0 

Q2IC=10000.      !IC  FOR  PRODUCER  STORAGE  !3.0 

147      CONTINUE 

Q3IC=50. 
C        THRESH=500.      'DOUBLING  THRESHOLD  FOR  CONSUMER 

IF  (IOFLAG.EQ.0)GOTO521 

CALL  GGON 

IF  (IPTR.EQ.0)GOTO151  13.0 

WRITE(3, 20171) 
151      CALL  GGINIT 

CALL  GGAXIS(0,0,767,479) 

CALL  GGERA 

CALL  GGBOX(7,0,0,767,479) 

WRITE(3,51) 
C 
C        CLEAR  MACROS  AND  DEFINE  ONE  TO  DRAW  BOXES 

51  FORMAT{ '+','9.  I?:A  P(+0,+01W(S1)V[,+24]V[+24,lV[,-24]Vt-24,] 
+   W(SO)  9;') 

WRITE    (3,511) 
511  FORMAK '+' ,'§:B   T(A1)    P[+0,+0]  V[  ,+24]  V[+24  ,]  V[  ,-24]  V[-24  ,] 

+     W(S0)T(A0)    a,-') 
WRITE(3,551) 

551  FORMAT( '+L(A1) •/ 
+'+L"7"FFFFFFFFFFFFFFFFPFPF; '/ 
+'+L"6"AA55AA55AA55AA55AA55; '/ 
+'+L"5" 92492492492492492492; '/ 
+  '+L"3"84210842108421084210;'  ) 

WRITB(3,552) 

552  FORMAT( ' +L"4"8844221 18844221 18844; ■/ 
+'+L"2"42009 10024009 1004200; '/ 
+'+L"1"200 0084002100 00420 00;'/ 
+'+L"0" 00002000000200 002000; '/ 

+  '+L"B"O00000OO000O00O00000;  '  ) 
DO  52  1=0,7 

CALL  GGPLT(I,735,(I+1)»24-16,0) 
CHAR=I+48 

52  WRITE(3,398)CHAR 

CALL  GGB0X(7, 725, 0,767, 248) 

CALL  GGBOX(7,0,0,767,248) 

CALL  GGB0X(7, 575, 0,725, 248) 
521      CONTINUE 
C 
C 
C 
C 


234 


SET  OP  SURFACE  OF  FORCING  ENERGY 


DATA  TEXT/8 0*0/ 


DATA  E  /O 

0,0.8,0.8 

0,0.8 

0,0.8 

0,0.8 

0,0.8 

0,0.8 

0,0.8 

0,0.8 

0,0.8 

0,0.8 

0,0,0 


0,0,0,0,0,0,0,0,0,0,0, 

0.8,0.8,0.8,0.8,0.8,0.8, 

1,1,1,1,1,1,0.8,0, 

1,1,1,1,1,1,0.8,0, 

1.2,1.2,1.2,1.2,1,1,0.8,0, 

1.2,1.2,1.2,1.2,1,1,0.8,0, 

1.4,1.4,1.4,1.4,1,1,0.8,0, 

1.2,1.2,1.2,1.2,1,1,0.3,0, 

1,1,1,1,1,1,0.3,0, 

1,1,1,1,1,1 ,0.8,0, 

0.8,0.8,0.8,0.8,0.3,0.8,0.8,0.8,0, 

0,0,0,0,0,0,0/ 


0.8,0.8,0, 
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X1(1)=5 

Y1(1)=5 

Q3( 1)=Q3IC 

ICON(5,5)=1 

IXYZ( 1)=IXY(5,5) 

K1=.417E-5 

K2=.5 

K3=.05 

K4=.45 

K5=.5E-4 

K6=.45E-3 

K7=.2E-6 

K8=.18E-5 

K9=.2E-5 

K10=.417E-5 

K11=.5E-3 

K12=.05 

K13=.7833E-6 

DO  200  1=2,11 

DO  200  J=2,11 

Q1{I,J)=Q1IC 

Q2(I,J)=Q2IC 

Q4(I,J)=Q4IC 

E(I,J)=E(I,J)«100. 

CONTINUE 

DO  201  IK=1,12 

Q4( 1 ,IK)=Q40IC 

Q4( 12,IK)-Q40IC 

Q4(IX,1)=Q40IC 

Q4(IK,12)-Q40IC 

CONTINUE 


!SET  PREDATOR  CLOSE  TO  CENTER 


1  FIRST  PREDATOR  LOCATION 
1 CODED  LOCATION 


!DK 
!DK 
!DK 
!DK 
!DK 
!DK 


CHANGE  ENERGY  LEVEL? 

IF( lETYP.EQ. 1 )GOT0221 
DO  220  1=2,11 
DO  220  J=2,11 
E(I,J)=(XMEAN*RAN(-1)+.5)*100. 


! SCALE  RANDOM  FUNCTION 
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220 
221 


270 
C 

C 

c>>>>>> 

c 

300 


IF  (IETYP.EQ.3)GOTO220 

E(I,J)-XMEAN*100. 

CONTINUE 

DO  250  1=2,11 

DO  250  J=2,11 

ETOT=ETOT+B(I,J) 

CONTINUE 

IF  (IETYP.EQ.2)GOTO300 

SF=ETOT/( 100. •100.) 

ETOT=0.0 

DO  270  1=2,11 

DO  270  J=2,11 

E(I,J)=E(I,J)*XMEAN/SF 

ETOT=ETOT+E(I,J) 

CONTINUE 


! CHANGED  3  TO  2  IN  3.0 


LOOP  START  <<<<<< 


CONTINUE 

IT=T 

EUSED=0.0 

SPTEMP=0. 

PTEMP=0 . 

DO  400  1=2,11 

DO  400  J=2,11 

RATE  EQUATIONS 

17=0 

ItIUM=0 

XEQ=0.0 


!  LOOP  START 

!GET  INTEGER  VALUE  OP  TIME  FOR  MOD  FUNCTION 

'ENERGY  USED  PER  TIME  I.E.  POWER 


IF    (ICON(I,J) .NE.0)XEQ=1.0 

J0=E(I,J) 

JR=JO/( 1+K13*Q1(I,J)*Q4(I,J)) 

R1=DT*K1*Q1(I,J)*Q4(I,J)*JR 

R2=DT*K2*Q1(I,J) 

R3=DT*K3*Q1(I,J) 

R4=DT*K4*Q1( I, J) 

R5=DT*K5*g2(I,J) 

R6=DT*K6*Q2(I,J) 

R10=DT*K10*Q1(I,J)*Q4(I,J)*JR 

R11=DT*K11*Q2(I,J) 

EUSED=EUSED+( JO-JR) *DT 

LEVEL    EQUATIONS         .... 


PTEMP=PTEMP+R1 


! PRIMARY    PRODUCTION 


Q1(I,J)=Q1(I,J)+R1-R2 

IF(Q1(I,J).LT.0.0)Q1(I,J)=0.0 

Q2(I,J)=Q2(I,J)+R3-R11 

Q4( I,J)=Q4( I,J)+R4+R6-R10 

Q4(I,J)=Q4(I,J)+R5*( 1-XEQ) 

ADD  LINEAR  FLOW  TO  Q4  IF  Q3  NOT  THERE 
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c 

C...         CONSUMER  CHECKING    ROUTINE 

C 

IP(XEQ.EQ.0)GOTO350  ! SKIP  IF  NO  CONSOMER  PRESENT 

IXYLOC=IXY(I,J)  ! GET  CODED  LOCATION 

INUM=1  lAT  LEAST  ONE  CONSUMER  PRESENT 

DO  217  I7=1,T2  !GET  CONSUMER  NUMBER 

IF  (IXYZ(I7).NE.IXYL0C)  GOT0217  IWRONG  CONSUMER  GOTO  217 

C...     RATE  EQUATIONS  FOR  CONSUMERS 

C 

XDT=DT 

C...     ...IF  RATIO  OF  Q2/Q3  IS  TOO  LOW  THEN  ITERATE  MORE  SLOWLY 

IF(Q2(I,J)/Q3(I7).LT.5.0)XDT=.01 
DO  650  DDT=XDT , DT , XDT 
R7=XDT*K7*Q2( I , J) *23 ( 17 ) *Q3( 17 ) *XEQ 
R8=XDT*K8*Q2 ( I , J) 'QS ( 17 ) *Q3 ( 17 ) *XEQ 
R9=XDT*K9*Q2 ( I , J) 'QS ( 17 ) *Q3 ( 17 ) *XEQ 
R12=XDT*K12*Q3(I7)*XEQ 

SPTEMP=SPTEMP+RS*(XDT/DT)+R7  1 SPTEMP  =  CONSUMP 

C 

C...  LEVEL    EQUATIONS    FOR  CONSUMERS 

Q3 ( 17 ) =Q3 ( 17 ) +R5/ ( ICON( I , J) ) * ( XDT/DT) +R7-R1 2 
Q2(I,J)=Q2(I,J)-R9  I  UPDATE   PREY   CONSUMED 

Q4(I,J)=Q4(I,J)+R8+R12  'UPDATE   NUTRIENTS 

650  CONTINUE 

C 

C 

C 

IF((T-XTIME)  .LT.1.)GOT0217  !  SKIP  MOVEMENT  IF  NOT  WHOLE  DT 
XTIME=T 
C        CHECK  PRESENT  POSITION  FOR  VALUE  OF  Q2 
QMAX=0 

IF(Q2(I,J).LT.PTHRSH)G0T0  1457 
XTEMP=X1(I7) 
YTEMP=Y1 ( 17 ) 
GOTO600 
C   IF  Q2  IS  STILL  CONSUMABLE  DON'T  MOVE,  JUST  EAT  SOME  MORE 
1457     DO  600  I2=I-N,I+N 
DO  600  J2=J-N,J+N 
IF  (I2.LT.1)  GO  TO  600 
IF  (I2.GT.12)  GO  TO  600 
IF  ( J2.LT.1)  GO  TO  600 
IF  (J2.GT.12)  GO  TO  600 
IP(Q2(I2,J2).LT.QMAX)GOTO580 
QMAX=Q2(I2,J2) 
XTEMP=12 
YTEMP=J2 
580      CONTINUE 
600      CONTINUE 

IC0N(I,J)-IC0N(I,J)-1 
CI  REMOVE  CONSUMER  FROM  PRESENT  LOCATION 

ICON ( XTEMP , YTEMP ) =ICON ( XTEMP , YTEMP ) + 1 
CIMOVE  CONSUMER  TO  NEW  LOCATION 
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CI CODE 
C 

C... 
C 


2000 

C 

C 


c 
c 

2001 
C 


IXYZ ( 17 ) =IX Y( XTEMP , YTEMP ) 
NEW  LOCATION 

CHECK  TO  SEE  IF  IT  IS  TIME  TO  REPRODUCE 


i IF  GREATER  THAN  REPSODUCTION 
'INCREASE  NUMBER  OP  CONSUMERS 
! ALLOW  NO  MORE  THAN  100 


IF(Q3(I7) .LT.THRESH)GOTO2000 

T1=T1+1 

IF(T1.GT.100)GOTO2000 

S3(I7)=Q3(I7)/2. 

Q3(T1)=Q3(I7) 

X1(T1)=XTEMP 

Y1(T1)=YTEMP 

IXYZ(T1)=IXy(XTEMP,YTEMP) 

IC0N(XTEMP,YTEMP)=IC0N(XTEMP,YTEMP)  +  1 
CONTINUE 


X1(I7)=XTEMP 

Y1(I7)-YTEMP     I  REM  REME.MBER  WHERE  TO  START  NEXT  TIME 


350 
C 
C 
400 


CONTINUE 
CONTINUE 


CONTINUE 
GPP=GPP+PTEMP 


1 ACCUMULATE  TOTAL  GPP 


453 

C 

C 


T2=T1 

IF  (T2.GE.100)T2=100 


END  CONSUMER  LOOP  AND  DO  BOOKEEPING 

CNSUMP=CNSUMP+SPTEMP 
PROD=GPP  +  CNSUMP 
TOTPOW=TOTPOW+EUSED 

COUNT  UP  THE  CONSUMERS 

NPROD=100 

1P(T2.LT.100)NPROD=T2 

Q3TOT=0 . 

DO  453  I=1,NPR0D 

Q3TOT=Q3TOT+g3 ( I ) 

ICNT{X1(I),Y1(I))=ICNT(X1(I),Y1(I))+1 

CONTINUE 

COUNT   UP    PRODUCERS    AND   NUTRIENTS 

Q1TOT=0. 

Q4TOT=0 . 
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Q2TOT=0 . 

DO  4531  1X1=2,11 

DO  4531  1X2=2,11 

Q1T0T=Q1T0T+Q1(IX1 ,1X2) 

Q2T0T=Q2T0T+Q2(IX1 ,1X2) 

Q4TOT=Q4T0T+Q4 ( IX 1 , 1X2 ) 

4531  CONTINUE 
Q4TOUT=Q4( 1,1)*44. 

TOT=Q 1 TOT+Q2TOT+Q3TOT+Q4TOT 
IF(IOFLAG.EQ.0)GOTO5000 
C 

c 

C WRITE  TEMPORARY  INFORMATION  AND  PLOT  GPP,  POWER  (FUSED) 

C 
C 

IF  ((T-PTIMEI).LT.I .OjGOTO  5000 

PTIME1=T 

ENCODE (80, 2006, TEXT) T,T2,EUSED/( 100.*DT) , ETOT/1 00 ., VERS , TITLE 
S,ETYPE(IETYP) 
2006     FORMAT( 1X,'T=' ,F6.2,'  C0NS=',I3,'  POW  USED= ' , F6 . 2 , 
S   '  AVAIL  P0t'J=  ',F6.2,'  VER:  '  ,P5  .2,  IX,  10A1  , 1 A4) 

CALL  GGTEXT(7, 0,475, TEXT, 1,0) 

IPT=PTEMP/(2.*DT*1000. ) 

CALL  GGPLT(4,IT,IPT+250,1)  ! GREEN  =  GPP 

IP0WER^( (EUSED/( 100.*XMEAN»DT) )-80.)*5.  lOUTPUT  80  TO  100 

CALL  GGPLT(2,IT,IPOWER+250,1)  ! RED  =  POWER 

C 

ENCODE (80, 4532, TEXT) Q1T0T,Q2T0T,Q4T0T,Q3T0T, TOT 

4532  FORMAT( 1X, 'Q1=  ',F10.2,'  Q2=  ',P10.2,'  Q4=  ',F10.2, 
S   '  Q3T0T=  ',F10.2,'  TOT=  ',P10.2) 

CALL  GGTEXT(6, 0,460, TEXT, 1 ,0) 

IYT=Q2TOT/20000. 

CALL   GGPLT(3,IT,250+IYT,1)  IMAGENTA  =    PRODUCERS 

IYT=Q3TOT/200 .+250 . 

CALL  GGPLT(7,IT,IYT,1)  ! WHITE   =CONSUMERS 


c 
c 

DRAW   PRODUCERS 

c 

c 

IF( (T-PTIME) .LT.TINT)GOTO4500 

PTIME=T 

DO    4050    1=2,11 

DO   4050    J=2,11 

I 
1 

1 

! 

C0L0R=Q2( I, J) 7(2*1000.) 

!3.0    CHANGED    4*1000    TO    2*1000 

IF(COLOR.GT.7)COLOR=7 

! 

1X9=1*24-40 

! 

JY=J*24-44 

1 

CHAR=48+C0L0R 

J 

IF(CHAR.GT.57)CHAR=57 

! 

CALL   GGPLT(COLOR,IX9,JY,0) 

(POINT   TO   LOWER  LEFT 

HRITE(3,39a)CHAR 

! 

398 

PORMAT( •+• , •T(A1)W(S' ,1H' ,A1,1H' 

,')    @B')                       ! 

399 

FORMAT( '       @A' ) 

(DRAW   BOX    (MACRO) 

4050 

CONTINUE 

! 
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3991 

C 

c 

c 

4500 
C 

c»>>> 

c- 
c 


CHAR='B' 

WRITE  (3, 3985) CHAR 

CALL  GGPLT(0, 600, 10,0) 

CALL  GGBOX(0,600,10,700,230) 

CALL  GGFILL(O) 

CALL  GGBOX(7,600,10,700,110) 

CALL  S0RT1(Q2, IX, 144,400.) 

CALL  GGPLT(7,601,IX( 1)+10,1) 

DO  3991  IT1=2,100 

CALL  GGVEC(3,600+IT1,IX(IT1)+10) 

CONTINUE 


CONTINUE 
>  PLOT  CONSUMERS  «<<<< 


ICOLOR=0 

IF(  (T-PTIMEC)  .LT.TINTOGOTO   5000 

PTIMEC=T 

CHAR='B' 

WRITE    (3,3985)CHAR 

FORMAT( '+' , 'T(A1)W(S'  ,1H'  ,A1 ,1H' , ' ) ' ) 

CALL  GGBOX(ICOLOR,284,4,560,244) 

CALL   GGFILL(O) 

CALL  GGPLTI 7,0,0,1) 

CALL  GGVEC(7, 767,0) 

DO    2005    11=1, T2 

ICOLOR=Q3(I1)/50. 

IF(ICOLOR.GT.7)ICOLOR=7 

CHAR=ICOLOR+48 

CALL  GGPLT(ICOLOR,X1(I1)*24-4+284,Y1(I1)*24-44,0) 

WRITE(3,398)CHAR 

CONTINUE 

CALL  GGB0X(7, 600, 130,700,230) 

CALL  SORT1(Q3,IX,T2,40.) 

CALL  GGPLT(7,601,IX(1)+130,1) 

DO  2017  IT1=2,T2 

CALL  GGVEC(7,600+IT1 ,IX(IT1 )+130) 

CONTINUE 
IF(IPTR.EQ.0)GOTO5000 
WRITE  (3,20171) 
FORMAT( '+S(H) ') 


2005 


2017 


20171 
C 

C 

C 

C 

5000 

C 

C 

C 


SEE  IP  ITS  TIME  TO  QUIT 
CONTINUE 

DIFFUSION 

QXT=0 . 0 


!DK 
!DK 
IDK 
!DK 
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DO  5002  1=2,11  IDK 

DO  5002  J=2,11  IDK 

DO  5001  l'P=I-1,I+1  IDK 

DO  5001  JT=J-1,J+1  IDK 

QXT=QXT+DK*(Q4(IT,JT)-Q4(I,J) ) ♦DT  !DK 

5001  CONTINUE  IDK 
Q4T(I,J)=Q4(I,J)+gXT  IDK 
QXT=0.0  IDK 

5002  CONTINUE  IDK 
DO  5003  1=2,11  IDK 
DO  5003  J=2,11  IDK 
Q4(I,J)=Q4T(I,J)                                          IDK 

5003  CONTINUE  IDK 
T=T+DT 

IP(T.LT.TTIME)GOTO300 
C>>>>>>  END  OF  MAIN  LOOP  <<<<<< 
C 

CALL  GGOFF 

DO  439  1=1,12 

DO  439  J=1,12 

ICNT( I , J) =ICNT( I , J) *DT 

439  CONTINUE 

CALL  ASSIGN(4, 'SURF4' ) 

WRITEC 4, 440) VERS, TITLE, BUF1,BUF2 

440  FORMAT( ' 1' , 'SURFACE  MODEL  VERSION  NO.  ' ,F6 .2 , IX, 10A1 , 
&   1X,9A1,1X,3A1) 

WRITE( 4,454) ETOT,ETyPE( lETYP) , PROD, TOTPOW, TOTPOH/ ( TTIME* 100 . 
&  GPP,CNSUMP,N,ISSUC,DK 

454  F0RMAT(1X,'  INPUT  ENERGY  TOTAL=  ',F10.2,'  ENERGY  TYPE  ',A4/ 
&   IX,'  TOTAL  PRODUCTION  =',G15.6/ 

&   IX,'  TOTAL  POWER  USED  =',G15.5,'  AVE  POWER/CELL  =  ',G15.5/ 

&   IX,'  GPP=  ',G15.6,'  TOTAL  CONSUMPTION=  ',G15.6/ 

S   1X,'  SEARCH  LENGTH     =',I3,'  STARTING  CONDITION  =  ',12/ 

5  IX,'  DIFFUSION  COEFFICIENT  =  'F7.5) 

WRITE( 4,455 )TTIME,DT,Q4 TOT/ 1000. ,Q4TOUT/1000 . , 

6  (Q4TOUT+Q4TOT)/1000. 

455  F0RMAT(1X,'    FOR    ' ,F10 .0 ,' ITERATIONS      DT=    ',F7.3/ 
S      IX,'  TOTAL   NUTRIENTS    (KG)       =    ',F10.2/ 
S       11X,'Q4   OUTER  TOTAL    (KG)    =    'F10.2/ 

S       11X, 'TOTAL   INNER  AND   OUTER    (KG)    =    ',F10.2/ 

5  1X,10X,'    NUTRIENT   MATRIX    Q4(I,J)'/) 

WRITE( 4,456) ( (Q4( I , J)/ 1000 . ,1=1 , 12 ) , J= 12 , 1 ,- 1 ) 

456  FORMAT( 1X,12F7.2) 

WRITE (4, 457 )PTHRSH, THRESH, Q2T0T/ 1000. 

457  FORMAT (// IX, 'VALUES    F0RQ2(I,J)    PRODUCERS'/ 

6  IX, 'PRODUCER  THRESHOLD  FOR  CONSUMER  MOVING=  ',F10.2/ 
&  IX, 'CONSUMER  THRESHOLD  FOR  DIVIDING  INTO  =  ',F10.2/ 
S       1X,'  TOTAL   PRODUCERS    (KG)       =    ',F10.2) 

WRITE( 4,4561 )( (Q2( I , J)/ 1000 . ,1=2 , 1 1 ) , J=1 1 , 2 ,-1 ) 
4561  FORMAT (8X,10F7. 2) 

WRITE( 4,4581) 

4581  FORMAT(//'    CONSUMER  VISITATION   MATRIX    '/) 
WRITE( 4,4582 )( (ICNT(I,J) ,1=1 , 12 ) , J=12 , 1 ,-1 ) 

4582  FORMAT( 12( 1X,I6) ) 
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WRITE(  4,4584) 

4584  FORMATC//'       FINAL   CONS'JMER  DISTRIBUTION' ) 
WRITE(4,45a5)( (ICON(I,J) ,1=1 , 12) , J=12 , 1 ,- 1 ) 

4585  FORMAT( 12( 1X,I6) ) 
WRITE(4,4571)Q3TOT 

4571  FORMAT(// IX, 'VALUES    FOR  CONSUMERS   TOTAL   CONSUMERS   =',G15.6) 

MRITE(4,458)((I,Q3(I),X1(I),Y1(I),IXYZ(I)),I=1,NPRDD) 
458  FORMAT ( IX, 14,'    Q3=    ',F8.2,'    X=',I2,'       Y  =" ,12 , IX ,14 ) 

CALL   CLOSE( 1 ) 

END 
C 
C 
C 

SUBROUTINE   S0RT1 ( X,IX ,N,SF) 

DIMENSION   X( 1) ,IX( 1) 

DO   20    1=  1  ,N 

IX(I)=X(I)/SF 
20  CONTINUE 

DO   40    1=1, N 

DO    40    J=I,N 

IF(IX( J).LT.IX(I))GOTO40 
ITEMP=IX(I) 
IX( I)=IX( J) 
IX( J)=ITEMP 
40  CONTINUE 

RETURN 

END 
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PROGRAM  GRAPH2 
C        vers  RGRF 
C 

C        WRITTEN  BY  JOHN  RICHARDSON 
C 

C        CALL  TO  REGLIN  ADDED  6/15/83 
C 

C        MAXIMUM  NUMBER  OF  DATA  POINTS  SET  TO  250 
C 

c 

C        MODIFIED  4/11/83  FOR  RGL  LIBRARY 
C 

BYTE  XTEXT(80) ,YTEXT(80) ,TITLE(80) ,ESC 

BYTE  PNAME( 16) ,IFNAM( 16) 

LOGICAL  I RPL AG, SMOOTH, SHADE 

DIMENSION  X(250) ,Y(250) ,Y1(250) 

DATA  TITLE/80*0/ 

DATA  XTEXT/80*0/ 

DATA   YTEXT/80*0/ 

DATA   FNAME/16*0/ 

DATA   IFNAM/16*0/ 

ESC=27 

TYPE    10 

10  FORMATC    GRAPHING    PROGRAM   FOR  GIGI' 

&/'  COMPLIMENTS   OF   JOHN    RICHARDSON'/) 

WRITE(5,105)ESC,ESC 
105  F0RMAT(2X,A1, 'PrTMV  ,A1, '«' ) 

TYPE    1  1 1 
111  FORMAT ( '    PROGRAM   REQUIRES   THE  TT:    BUFFER   BE   SET   TO   NOWRAP'// 

S'    SET   /NOWRAP=TI : ' // 

S'    IF   THIS    IS   NOT   DONE   PLEASE   EXIT   PROGRAM   AND   CORRECT   THIS'//) 
TYPE    15 

15  FORMATC       IS    DATA   IN    A   DATA   FILE?    (1=YES,    0=NO,    -1=EXIT)     '$) 
ACCEPT    16,IANS1 

16  F0RMAT(I2) 

IF    (IANS1 .LT.O) STOP 'MAKE   CHANGES    AND    RERUN' 
IF    (IANS1.EQ.0)    GOTO    11 
TYPE    161 

161  FORMAT( '    FILE   NAME    FOR  DATA:     '$) 
ACCEPT    162,     (FNAME(I) ,1=1 ,16) 

162  FORMAT! 16A1) 
0PEN(UNIT=1,NAME=FNAME,TYPE="0LD'  ,FORM=' FORMATTED'  ) 
READ    ( 1 ,1621) TITLE 

1621  FORMAT( 1x,80A1) 
HRITE( 5, 1621) TITLE 
READ(1,1622)NPAIRS 

1622  FORMAT( 1X,I3) 

DO  1630  I=1,NPAIRS 

READ(  1,*)X(I)  ,Y(I) 
C1625  FORMAT(2G15.6) 

1630  CONTINUE 

N=NPAIRS 

GOTO  51 

11  TYPE  20 
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20       FORMAT ( '   HOH  MANY  PAIRS  OF  POINTS  TO  PLOT  '$) 

ACCEPT  30,N 
30       FORMAT (13) 

IF(I.GT.500)GOTO11 

TYPE  35 

35  FORMAT ( IX, 'DESCRIPTION   OF   DATA    (UP   TO   SO    CHARACTERS    '/) 
ACCEPT   36,(TITLE(K) ,K=1 ,80) 

36  PORMAT(80A1) 
DO    50    1=1 ,N 

59  TYPE   60,1 

50       FORMAT ( '   X  AND  Y  VALUES  FOR  POINT  ',13, 
S '  SEPARATED  BY  COMMAS  [ R]  ' ) 
READ  (5,*,ERR=9911)X(I) ,Y(I) 

50  CONTINUE 
CLOSE (0NIT=1) 

WRITE( 5,601) 

601  P0RMAT(//1X,'D0   YOU  WANT  TO   SAVE  DATA    ( 1-YES,    0-NO)'$) 
READ(5,6D2)ISAVE 

602  FORMATdD 

IF    (ISAVE.NE.1)G0T051 
HRITE(  5,603) 

603  FORMAT ( '       WHAT    IS   THE   FILE   NAME    FOR  THE   DATA    '$) 
READ(5,604)IFNAM 

604  P0RMAT(16A1) 

CALL    ASSIGN(2,IFNAM) 

HRITE(2,606)TITLE 
606  PORMAT( 1X,80A1) 

WRITE(2,608)N 
608  FORMAT( 1X,I3) 

DO    511    1=1, N 

WRITE(2,611)X(I)  ,Y(  I) 
611  FORMAT(1X,G15.6, '     ,     ',G15.6) 

511  CONTINUE 

CL0SE(UNIT=2) 

51  XMAX=a. 
YMAX=0. 
A=0. 
3=0. 
R2=0. 
CEE=0. 
TYPE   5001 

5001  FORMAT ( '  DO  YOU  WANT  TO  RUN  REGRESSION  ON  DATA?  (1-YES,  0-NO) ' ) 
ACCEPT  5002, IRGS 

5002  FORMATdD 

IF    ( IRGS. NE.O) CALL    REGLIN( N,X, Y, A,B ,R2 ,CEE) 
TYPE   501 

501  FORMAT( '       WHAT    IT   THE   X-    AXIS   DESCRIPTION') 
ACCEPT    502,XTEXT 

502  FORMAT(80A1) 

CALL   STRIP(XTEXT,80) 
TYPE    503 

503  FORMAT('    WHAT    IS   THE   Y-AXIS   DESCRIPTION') 
ACCEPT    504,YTEXT 

504  PORMAT(80A1) 
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CALL   STRIP(5fTEXT,80) 

TYPE   70 
70  FORMAT!'    WANT   TO    INPUT   MINIMUMS    AND   MAXIMUMS    (1-yes,    0-no)'$) 

ACCEPT    72,I!4IN 
72  FORMAT(II) 

IF    (IMIN.EQ.0)GOTO85 

TYPE   74 
74  FORMAT(  1X,"WHAT    ARE   XMIN   AND   XMAX    [  R]     'S) 

ACCEPT    *, XMIN, XMAX 

TYPE   76 
76  FORMAT( IX.'WHAT    ARE   YMIN    AND    YMAX    [ R]     '$) 

ACCEPT    *, YMIN, YMAX 
85  CONTINUE 

TYPE   89 
89  F0RMAT(1X,'    LINE  TYPE    (0-9)     '$) 

ACCEPT    891 ,ILIN 

891  FORMAT(II) 
8910           TYPE   892 

892  F0RMAT( IX, 'VALUE    FOR  DATA  MARKER    (0-9,    -1    TO    SEE   LIST)     '$) 
ACCEPT    893,IMAR:< 

IF(IMARK.GE.0)GOTO8930 
WRITE( 5,8921) 
8921  FORMAT(/'     0-    POINT' 


& 

/'     1-    SQUARE' 

& 

/'    2-    OCTAGEN' 

& 

/'     3-   TRIANGLE' 

S 

/'    4-   CROSS' 

& 

/'     5-   X' 

s 

/'    6-    Y' 

s 

/'    7-    DIAMOND' 

& 

/'    8-   ARROWHEAD' 

& 

/'    9-    HOURGLASS' 

& 

/'     10-POINT    IN   A  CIRCLE') 
GOT089 1 0 

8930 

CONTINUE 

393 

FORMAT ( 15) 

IF    (1MIN.EQ.1)IROUND=0 
IF    (IMIN.EQ.1)GOT09910 
TYPE   90 

90 

FORMATC       ROUND   MAX    AND   MIN  VALUES?    (1-YES,    0-NO) 
ACCEPT   99,IR0UND 

99 

FORMAT(II) 

IF(  IROrjND  .GT .  1  )  G0T05 1 

1F(IROUND.LT.O)GOT05  1 

9910 

IRFLAG=. FALSE. 

IFdROUND.EQ.I  )  IRFLAG=  .TRUE. 
WRITE    (5,9901) 

9901 

FORMATC       CURVEFIT   THE   DATA  LINE    (1-YES;    0-NO) '5) 
READ(5,9902)ISM 

9902 

FOPMAT(II) 

SMOOTH=. FALSE. 

IF(  ISM. EQ.1)SM00TH=. TRUE. 

WRITE(5,991)ESC 

991 

FORMAT( '+' 1A1 ,' tH' )                                    ! SEND   CURSOR  HOME 

'$) 
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SHAD B= . FALSE . 

CALL  INITGR(S) 

CALL  CLRSCR 

CALL  CLRTXT 

CALL  SCOLOR( ' GRAYO ' , 0 ) 

CALL  SC0LOR( ' GRAY 1 ' , 1 ) 

CALL  SCOLOR( ' GRAY2  '  , 2 ) 

CALL  SCOLOR( ' GRAY3 ' , 3 ) 

WRITE (5, 991) ESC 

CALL  DPAPERCLIN' ,10,2,'LIN' ,10,2,"GRAY3' ) 

IF  (IMIN.EQ.0)GOTO1211 

CALL  LNAXIS{ 'YL' ,yTEXT,YMIN,YMAX,IRPLAG) 

CALL  LNAXIS( 'XB' ,XTEXT,XMIN,XMAX,IRFLAG) 

GOT01212 

1211  CALL  LNAXISC 'YL' ,YTEXT,,,IRPLAG) 
CALL  LNAXIS( 'XB' ,XTEXT, , ,IRFLAG) 
WRITE (5,991) ESC 

1212  CALL  PDATA(N,X,Y,'L','GRAY2',IMARK,ILIN, SMOOTH, SHADE, 0.0) 
IP(IRGS.Eg.0)GOTO1234 

DO  1277  11=1, N 
Y1(I1)=B*X(I1)+A 

1277  CONTINUE 

CALL   PDATA(N,X,Y1,'L' , 'GRAYS' ,0 , 1 , .FALSE. , .FALSE. ,0 .0) 

TYPE    121, ESC 

WRITE    (5,1278)B,A,R2,CEE 

1278  F0RMAT(/1X,'    Y=    ',G12.4,'*X+    ',G12.4,'     :R92    =    ',G12.4, 
+'STD    ERR  =    ' ,G12.4) 

1234  TYPE    121, ESC 

121  F0R!1AT('+',1A1,'[H   0-QUIT;     1-REPLOT;    2-SCREENDUMP' ,  $ ) 
ACCEPT    122,    lANS 

122  F0RMAT(I2) 
IF(IANS.EQ.1)GOT0   51 

IF    (IANS.NE.2)G0T0   2550 

WRITE(5,1221)ESC 
1221  PORMAT( '+' ,1A1, ' [H' ,80X) 

CALL   CPYSCR 

GOT01234 
2550  STOP    'END' 

9911  WRITE( 5,9912) 

9912  FORMATC       ERROR   IN    ENTRY   PLEASE   RE-ENTER') 
GOTO   59 

END 


SUBRODTINE   STRIP    (TEXT,N) 

BYTE   TEXT( 1) 

DO   20    I=N,1,-1 

IF    (TEXT(I).EQ.32.0R.TEXT(I) .EQ.0)GOTO20 

TEXT(I+1)=0 

RETURN 
CONTINUE 
RETURN 
END 
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C 

c 
c 

SUBRODTINE  REGL1N(N,X , Y,A,B,R2 ,CEE) 
C 

C        BASED  ON  PROGRAM  IN  'COMMON  BASIC  PROGRAMS'  BY 
C        LON  POOLE  AND  MARY  BORCHERS  P.  145 
C 


DIMENSION  X( 1) ,Y( 1) 
REAL  J,K,L,M 


J=0.0 
K=0.0 
L=0.0 
M=0.0 
A=0.0 
B=0.0 
R2=0.0 
CEE=0.0 
DO  100  1=1, N 
J=J+X(I) 
K=K+Y(  I) 
L=L+X(I)*X(I) 
M=M+Y(I)*Y(I) 
R2=R2+X(I)*y(I) 
100      CONTINUE 
XN=N 

B=(XN*R2-K*J)/(XN*L-J*J) 
A=(K-a*J)/XN 
J=B*(R2-J*K/XN) 
M=M-(K**2)/XN 
K=M-J 
R2=J/M 

CEE=SQRT(K/(XN-2) ) 
RETURN 
END 
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